The State of the World’s Beaches

Coastal zones constitute one of the most heavily populated and developed land zones in the world. Despite the utility and economic benefits that coasts provide, there is no reliable global-scale assessment of historical shoreline change trends. Here, via the use of freely available optical satellite images captured since 1984, in conjunction with sophisticated image interrogation and analysis methods, we present a global-scale assessment of the occurrence of sandy beaches and rates of shoreline change therein. Applying pixel-based supervised classification, we found that 31% of the world’s ice-free shoreline are sandy. The application of an automated shoreline detection method to the sandy shorelines thus identified resulted in a global dataset of shoreline change rates for the 33 year period 1984–2016. Analysis of the satellite derived shoreline data indicates that 24% of the world’s sandy beaches are eroding at rates exceeding 0.5 m/yr, while 28% are accreting and 48% are stable. The majority of the sandy shorelines in marine protected areas are eroding, raising cause for serious concern.

server side of the platform reduces image processing time to only several minutes per image 13 and enables efficient validation of the automatically detected shorelines at multiple sites where ground-truth field data are available.
To enable global mapping of sandy shorelines it is first necessary to identify sandy beaches and then determine shoreline positions in every image in the GEE platform. The spatio-temporal scales associate with this study (i.e. global scale, 33 year analysis) and the large amount of satellite images that therefore need to be analysed necessitates the use of robust automated image analysis techniques. Machine learning 14 and image processing 15 techniques that lend themselves to such automated analyses are readily available. However, to be able to use satellite derived shoreline positions for real-world applications such as reliably estimating trends and structural damage to infrastructure, a horizontal resolution of at least 10-20 m is required. For example, shoreline change rates above 0.5 m/yr over a long period are typically employed to flag a coastal area as one experiencing chronic (=long term i.e. decades to centuries 7 ) erosion or accretion. Over a period of 30 years that would mean a total displacement of just 15 m. Previous studies have evaluated the positional accuracy of satellite derived shorelines (SDS) based on single images [15][16][17][18][19] to range between 1.6 and 10 m. It should be noted that these studies suffered from limitations such as the number of images used, the quality of the in-situ data used for validation or the magnitude of changes in observed shoreline position. Recently, Hagenaars et al. 20 presented a long-term, but local-scale satellite image analysis on shoreline trends, that overcomes all of the aforementioned limitations. They found the accuracy of the SDS derived from moving average composite images to be of subpixel precision (~half a pixel size, i.e., 15 m for Landsat and 5 m for Sentinel-2). The accuracy of <15 m, reported by Hagenaars et al. 20 for composite Landsat images, matches the required displacement of 15 m for reliable shoreline change classifications over the last 30 years. For that reason, we adopt the same approach in this study, yet at a global scale.
Here we present an up-to-date global-scale assessment of dynamics of sandy shorelines using a fully automated analysis of 33 years (1984-2016) of satellite images. First, we detect sandy beaches worldwide by applying a pixel-based supervised classification to a cloud-free high-resolution global composite image for 2016. A digital beach training dataset is provided to the classification software and validated for 50 locations worldwide that include both sandy and non-sandy beaches. Next, we apply a shoreline detection algorithm to cloud free global annual composite images using more than 1.9 million historical Landsat images. After a successful quantitative validation of this technique at multiple sites located in various geographical settings and environmental conditions, we derive shoreline change rates in m/yr at transects with an alongshore spacing of 500 m along the world's shoreline. The above mentioned methods are elaborated in the Methods section below while the complete validation is presented Supplementary Material (S2).
The main outcomes of our analysis include: (a) the global occurrence of sandy beaches, (b) rate of erosion/accretion at all sandy beaches in the world, (c) highlights of observed natural and human induced impacts on coastal erosion/accretion at selected locations, and (d) identification of global hot spots of coastal erosion/accretion.

Results
Global Occurrence of Sandy Shorelines. Coastal classifications have been widely employed in the field of geomorphology to characterise the diversity of coastal landforms and the contexts within which they emerge, but hitherto no single system of classification has been comprehensive in scope or coverage 21,22 . Criteria in these classifications typically include tectonic 23 and hydrodynamic controls, as well as the sedimentological response. Hydrodynamics controls considered include classifications of wave parameters 24 , tidal range 24,25 and a combination of both 26 . A ternary classification presented by Boyd et al. 27 , which considers the relative importance of fluvial inputs, wave energy, and tidal forcing provided a useful analysis of siliciclastic sedimentary coasts. The combination of tectonic and hydrodynamic controls led to the proposition of coastal morphogenetic classifications 28 , which are probably the most widely used classification schemes.
Sediment texture and composition 29 are additionally useful to classify and describe coastal sedimentary environments. However, previously reported values of the global occurrence of sandy shorelines vary between 10% and 75% (see Table 1). The methods used to arrive at these values remain, in most cases, unclear or qualitative (as also indicated in Table 1).
In our analysis, we applied supervised (human-guided) classification to global cloud-free satellite images (see Section 3.2) to identify sandy shorelines. One of the main reasons for our focus here on sandy beaches is that detecting shoreline dynamics for non-sandy shores like muddy coasts can be complex. Mild foreshore slopes, resulting in large horizontal tidal excursions, and high water content hampers correct shoreline detection. In the case of mangroves, seasonal growth cycles can impede correct shoreline detection. Moreover, it should be noted that as the reflectance signatures of sand and gravel beaches cannot be differentiated in the satellite imagery, as both materials originate from the same granular composites of finely divided rock, our references to sandy beaches herein also includes gravel beaches.
Our analysis showed that 31% of the ice-free world shoreline is sandy. The continent with the highest presence of sandy beaches is Africa (66%), while in Europe only 22% of the shoreline is sandy (see inserted table in Fig. 1). The percentage of sandy shorelines obtained from this analysis for USA and Australia compare well with the more recently reported regional scale values (see Table 1). The larger deviation in percentage found for Europe is significantly influenced by the smaller total length of shoreline used in the Eurosion 10 data base. It should be noted that the sandy beach classification also includes the gravel beaches in the world. The reflectance signatures of sand and gravel beaches cannot be differentiated in the satellite imagery as both materials originate from the same granular composites of finely divided rock.
The global latitudinal distribution of sandy shorelines shows a distinct relation with latitude and hence with climate; no relation is found with longitude. The relative occurrence of sandy shorelines increases in the subtropics and lower mid-latitudes (20°-40°) with maxima around the horse latitudes (near 30°S and 25°N; see Fig. 1). In contrast, they are relatively less common (<20%) in the humid tropics where mud and mangroves 30  Global sandy beach erosion. Worldwide beach erosion became apparent during the 1980s following the studies of the International Geographical Union working group on the Dynamics of Coastal Erosion (1972)(1973)(1974)(1975)(1976) and the Commission on the Coastal Environment (1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984). In these studies, two hundred participants representing 127 countries contributed to a survey which indicated that 70% (10%) of the world's sandy beaches experienced net erosion (accretion) while 20% were stable 32 . However, as these estimates were primarily a result of interviews, they are necessarily qualitative, at best. Furthermore, the estimates likely did not take into account changes occurring along undeveloped and uninhabited coasts due to the subjective methodology adopted.
The quantitative global distribution of sandy shorelines presented herein, for the first time, allows the derivation of objective and up to date global scale assessment of chronic shoreline changes (i.e. beach erosion/accretion). Beach erosion can occur at a range of timescales 33 . Individual storms will generally result in rapid short-term erosion, followed by short-term accretion, leading to negligible net change over time scales of a few weeks-months. If sediment deficiencies persist for long periods of time (e.g. due to longshore gradients in sediment transport, reduction of fluvial sediment supply to the coast), chronic erosion can result. The analysis presented here focusses on such chronic erosion and accretion. However, there are no common standards for the classification of rates of chronic beach change 34 which is generally quantified through some statistical treatment of erosion rates and/or volumetric losses (e.g. ref. 35 . The accuracy of the SDS data of ~0.5 pixel (see Section 1) and the study period of ~30 years allows for a classification of beach change rates with class boundaries of 0.5 m/yr. Hence, we adopted the chronic beach erosion classification scheme proposed by Esteves and Finkl 36 and extended it with a classification for extreme erosion resulting in the below scheme:  Table 1. Reported values of global and regional occurrences of sandy shorelines and percentages of chronic erosion and accretion.
Our assessment shows that 24% of the world's sandy beaches are persistently eroding at a rate exceeding 0.5 m/yr over the study period , while 27% are accreting (see Table 1). About 16% (18%) of sandy beaches are experiencing erosion (accretion) rates exceeding 1 m/yr.
Chronic erosion of beaches (<−0.5 m/yr) is shown across the globe with relatively low latitudinal variation (see Fig. 2). Generally, between 30% and 40% of sandy beaches per degree latitude are eroding with relatively high eroding values up to 50% just south of the equator associated with large-scale land losses adjacent to the Amazon River mouth.
More severe erosion rates are found at various locations across the globe. About 7% of the world's sandy beaches experience erosion rates classified as severe. Erosion rates exceed 5 m/yr along 4% of the sandy shoreline and are greater than 10 m/yr for 2% of the global sandy shoreline. On the other hand, about 8% of the world's sandy beaches experience significant accretion (>3 m/yr), while 6% (3%) are accreting more than 5 m/yr (10 m/yr).
Taking a continental perspective, Australia and Africa are the only continents for which net erosion (−0.20 m/yr and −0.07 m/yr respectively) is found, with all other continents showing net accretion. The continent with the largest accretion rate (1.27 m/yr; see table in Fig. 2) is Asia, likely due to the artificial development of the Chinese coast and large land reclamations in, for example, Singapore, Hong Kong, Bahrain and UAE. On a global scale, the world's beaches have accreted on average 0.33 m/yr over the past three decades, i.e. a total gain of 3,663 km 2 over this period.
Using the SDS data we then focussed on coastlines that are internationally recognised as nature protected areas by the World Database on Protected Areas (WDPA), which is the most comprehensive global database on terrestrial and marine protected areas, produced by UNEP-WCMC and IUCN 37 . Compared to the global average, a relatively high percentage of sandy shorelines in the WDPA-identified areas are experiencing erosion. Our analysis indicates that 32% of all marine protected shorelines are sandy of which 37% are eroding at a rate larger than 0.5 m/yr, while 32% are accreting.
Quantifying local scale erosion/accretion due to human interventions. No single explanation can easily account for the observed erosion/accretion trends along the global sandy shoreline, or for the acceleration of erosion/accretion on any particular beach 38 . However, analysis of local trends derived from the global scale shoreline assessment presented herein can help identify natural and human drivers of shoreline change. To illustrate this, we present two highlights of erosive behaviour and two of accretive behaviour. Another four highlights are presented in the Supplementary Material (S3). a) Sand mining and subsidence. The Mekong Delta in Vietnam, the third largest delta in the world, is increasingly affected by human activities and exposed to subsidence and coastal erosion. The large-scale shoreline erosion is attributed to excessive sand mining in the river and delta channels, and subsidence due to unregulated groundwater extraction 39 . Analysis of the SDS data (Fig. 3a) reveals slight erosion between 1984 and 1990, after which higher, but steady erosion rates are found. Erosion rates in the considered area typically range between 25-30 m/yr over the last three decades. Based on the strong linear trend, the SDS data may be used for projections of land loss and displacement strategies, as it is not expected that erosion rates will decrease in the near future unless mitigating measures are implemented.
b) Coastal structures. The harbour structures at Nouakchott, Mauritania, blocked the large unidirectional north-south longshore transport of sand since 1986, causing areas of beach erosion that has impacted the local social and urban developments. The shoreline evolution rates observed after the harbour construction are 10 times larger than the values that would have been observed in the natural state 40 . The harbour breakwaters induced severe erosion over a distance of more than 10 km in the downdrift zone where accretion was likely to occur in the absence of the harbour. The SDS data (Fig. 3b) shows erosion rates of 20 m/yr. c) Sand Nourishments. A large-scale bypass system became operational in 2001 at the Tweed River, New South Wales, Australia, to mitigate erosion of the beaches to the north of jetties constructed at the river entrance 41 . The bypass system pumps sand from south of the river mouth to three beach compartments located north of the river through buried pipelines. The SDS data (Fig. 3c) depicts a beach widening of ~250 m at Coolangatta Bay in the four years after the bypass system was commissioned. d) Interception of longshore drift by coastal structures. The construction of two training breakwaters at Praia da Barra near the Aveiro Lagoon, Portugal interrupted the high southward ambient alongshore transport estimated at about 1 million m 3 /yr 42 . This resulted in erosion at the south of the trained inlet affecting the shoreline over about 30 km downdrift, but also strong accretion updrift. The SDS data reveals the continuous and ongoing accretion of the northern beach at a rate of about 10 m/yr (Fig. 3d).

Global hot spots of erosive and accretive beaches.
Here we present the top eroding and accreting coastal stretches (i.e. hot spots) in the world ( Table 2). The largest erosive hot spot is just south of Freeport in Texas where a 17 km stretch the beach has eroded on average more than 15 m/yr over the last three decades. The world's longest coastal stretch suffering severe erosion is located farther to the east in Texas where we observed a 29 km stretch of sandy beach with a mean erosion rate of 5.3 m/yr. Interestingly, four of the seven largest hot spots are located in the USA, consistent with the widespread concern and reports of erosion in the USA 11,35,43,44 .
The largest accretive hot spot is in Namibia at a location where a mining company has built unprotected sandy bunds in the sea to facilitate the diamond prospecting. The area landward of the bunds is dried out to enable more convenient diamond prospecting. Naturally accreting beaches of lengths exceeding 20 km and change rates larger than 7 m/yr are found at a migrating barrier island (Schiermonnikoog, The Netherlands) and at locations where sand dunes migrate into the sea (Madagascar and Mauritania). It is noteworthy that four of the seven largest accretive hot spots are in fact human-induced.
Outlook. In the near future we foresee great potential for remote sensing techniques and big data analysis in operational monitoring of the World's coast and beaches. The global sandy shoreline change analysis presented herein is primarily based on Landsat imagery with 30 m resolution and a revisit time of 16 days. In recent years new satellites (Sentinel-2a,b) that will significantly enrich the satellite imagery data both in temporal (revisit time of a few days) and spatial resolution (<10 m) have been launched. At present, private institutions already provide satellite images at approx. 1 m resolution with a daily revisit and global coverage. We expect that this trend will demand more emphasis on big data statistics in the near future to closely and better monitor how the planet is changing.

Methods
The workflow applied in this study comprises three methods as discussed below and illustrated in Fig. 4.  previously reported values of 1 million km 3 , 1.16 million km 46 , and 1.47 million km 47 . In the future, we intend to merge the 500 m transect system with locally available grids and refine it where appropriate.

Global transect system.
Detection of sandy beaches. Sandy beaches are detected by applying a pixel-based supervised classification to a global Top of Atmosphere (TOA) reflectance percentile composite image for the year 2016 using all available Sentinel-2 images. To facilitate this, the world has been divided into boxes of 20 km × 20 km. Using the 2016 OSM shoreline, we only select the boxes that intersect with the 2016 shoreline, which results in about 24,000 boxes to be analysed. To train the supervised classifier, a beach area consisting entirely of sand is selected (at the Dutch Texel island) as well as training areas on land representing different types of land use. To select the most promising classification algorithm, the validation results were quantitatively compared to the sandy beach feature in OSM. From the four considered classification algorithms, the Classification and Regression Tree (CART) classifier resulted in the lowest omission error and the highest percentage of true positives (97%) using the beach features in a 100 km long section of sandy beaches along the Dutch coast.
Next, we apply the trained supervised classification method to all boxes to detect sandy beaches at global scale as the OSM beach feature is not available for the entire globe. A search area of 500 m land-and seaward of the 2016 OSM shoreline is defined, after which the supervised classification is conducted using GEE to automatically detect sandy beaches. The result is a series of polygons encapsulating all sandy beaches worldwide, including both quartz and carbonate sands, and gravel. More than 50 sand validation locations, randomly spread across the world, were selected independently from the training dataset. Validation through visual inspection resulted in 96% accuracy (see Supplementary Material S1).
Transects that intersect with a sandy polygon are classified as 'sand' and others as 'non-sand' . Transects for which no sand classification could be made due to the absence of a cloud-free Sentinel-2 image are labelled as 'undetermined sediment composition' . As this is applicable for 5.2% of all transects the percentage of sandy beaches is 31% ± 1.5%, assuming that the unknown areas behave similar to the global mean.

Dynamic shoreline detection.
To remove the effects of clouds, shadows, snow, and ice, we generate yearly top-of-the-atmosphere reflectance composites, which we then use to estimate an accurate surface water mask using dynamic thresholding method described in ref. 48 . Yearly composite images generated by the 15% reflectance percentiles per pixel were analysed to determine global shoreline positions, resulting in the removal of clouds and shadows. This approach is comparable to how Hansen 49 generates composite images. However, the use of an exact percentile value turns out to be more suitable than the interval mean averages used in that study. Analysis of the composite images significantly decreases the influence of the tidal stage on the detected shoreline positions and averages out seasonal variability in wave and beach characteristics. Nevertheless, at sites with persistent swell conditions the wave-induced foam due to wave breaking will introduce a seaward offset in detected shorelines. Fortuitously, however, this persistency ensures that the wave-induced offset is most likely also present in annual composites and shorelines of other years. Thus, the wave effects on detected shorelines are likely to be limited, especially where long-term shoreline change rates at such sites are concerned.
For validation purposes with long-term in-situ shoreline changes, an optimal averaging period of 192 days is applied; i.e. the first integer that is found when dividing the global revisiting time of the satellite sensor (16 days) by a semidiurnal tidal period (approx. 12 hrs). In case all satellite images in this averaging interval are cloud-free the average water level corresponds to mean sea level. The potential year-to-year random deviation from 'mean sea level' due to omitted satellite images is assumed to have a limited effect on the 33-year trend of shoreline change; this assumption will be verified as part of further research. Next, the resulting composite images are used to estimate the Normalized Difference Water Index (NDWI). The Canny edge detection filter is used to roughly estimate the position of the water-land transition, followed by the use of the Otsu thresholding method 50 on a buffer polygon around the water-land transition to identify the most probable threshold to classify water and land on the image. The detected water lines at the edge of the water mask are smoothed using a 1D Gaussian smoothing operation to obtain a gradual shoreline avoiding the pixel-induced staircase effect. A value of three gives the best results based on the four validation cases; meaning that it takes three cells on both sides during the 1D smoothing. The method may result in several shoreline vectors since lakes and small channels are detected. In this case, only the most seaward shoreline position is analysed. Other studies have applied global surface water change and occurrence detection 48,51 but they lack validation with in-situ measured shoreline changes. A number of studies have validated their methods with either cross-shore positions at one location 19 or over limited spatial scales 15 . Here we evaluate the validity of the shoreline detection method for four cases representing different types of beaches, sand, tidal and wave characteristics. Given the geographical spreading, we selected the following beaches with long-term shoreline monitoring programs: the Sand Engine (The Netherlands), Long Beach, WA (West Coast, USA), Narrabeen (Australia) and Hatteras Island (East Coast, USA). The latter case is presented below while the others are presented in the Supplementary Material (S2). For the same time period we collated 325 cloud-free satellite images and determined the shoreline position for this coastal stretch; the analysis took only 8 hours in total due to the computational power of the GEE platform. For each transect a linear regression was performed. The linear trends calculated from the SDS show good agreement with the observed shoreline change rates (see Fig. 5). The mean offset for all transects between observations and SDS is 2.0 m with a RMSE of 17 m.

Hatteras Island validation.
The Supplementary Material (S2) summarizes the error statistics for all four cases. Based on these validations, the shoreline detection method can be concluded to be capable in deriving long-term shoreline change rates for a variety of coastal settings. The average of the offsets over three validation sites is 2.3 m with a RMSE of 21 m.
Although the quantitative evaluation of the applied shoreline detection method with in-situ observations shows good capabilities, more verification is essential. Unfortunately, however, quantification of the influence of macro-tidal ranges, wave breaking and run up, beach slopes, etc. requires tidal, wave and beach characteristic information, which are generally not freely available.
Global change rates for sandy shorelines. For the global application presented here, we generated cloud-free annual-composites using the historical Landsat image archive. The automated shoreline detection method produces 33 annual global shorelines (1984-2016) with an alongshore resolution of 30 m. We then specified transects at a 500 m alongshore spacing, and determine the intersection point of each transect with the aforementioned annual shorelines, which provides a sequence of shoreline positions per transect. The shoreline change rate (m/yr) at each transect is then computed by applying linear regression to all shoreline positions at that location. Ideally, a SDS position is available for each transect annually. However, the availability of satellite images and cloud cover can limit the number of SDS positions. Encouragingly, however, 82% of all sandy transects consist of more than ten annual shoreline positions between 1984 and 2016. Nevertheless, to avoid unrealistic shoreline change rates we applied the following filters to all sandy transects: • Transects containing less than 5 (out of 33) SDS data points as well as transects with a temporal coverage shorter than 7 years are omitted from the analysis (9% of all transects).

•
Transects located beyond latitudes 60°N and 50°S (including Greenland and Antarctica) are omitted from the analysis due to possible ice coverage (9% of all transects).

•
In the linear regression, outliers are identified as SDS points deviating more than three times the standard deviation and hence not considered in the regression. If the remaining number of data points is smaller than 5 points, then the transect is omitted from the analysis. Applying these filters reduce the global data set to 81% of the original number of sandy transects. The linear regression method used to quantify long-term shoreline change rates performs well in capturing trends of chronic sandy shoreline change which is in line with findings of Crowell et al. 53 . However, multiple transects were characterised by unsteady changes in SDS positions for which other methods may be more appropriate. Ultimately, more than 60% of the 2.2 million transects show an uncertainty bandwidth of less than 50% of the linear trend rate, which can be considered as a proxy for the representativeness of the linear regression method.
The shoreline change rates, presented at an alongshore resolution of 500 m along the world's shoreline, will become publicly available and be accessible through the interactive website at: http://shorelinemonitor.deltares.nl.
Defining Hot Spots. In order to avoid localized hot spots, it was ensured that each eroding/accreting hot spot comprised at least 5 km of sandy shoreline where all considered transects showed either erosive or accretive change rates larger than 0.5 m/yr over the 33 year data set.
Two large-scale land reclamations appear in the top seven accretive beaches in the world. One reason is that those land reclamations consisted of bare sand in 2016, and hence are recognised as a wide sandy beach area by our methodology. The other reason is that the adjacent shorelines have advanced either due to the beach nourishment schemes or natural accumulation of sand in the shadow zones of these interventions.