Large regional variability in coastal erosion caused by ENSO

In the Pacic Basin, El Niño/Southern Oscillation (ENSO) is the dominant mode of interannual climate variability and drives substantial changes in oceanographic forcing, likely having a signicant impact on Pacic coastlines. Yet, how sandy coasts respond to these basin-scale changes has to date been limited to a few long-term beach monitoring sites, predominantly on developed coasts. Here we use 35 years of Landsat imagery to map shoreline variability around the Pacic Rim (72,000 beach transects) and identify coherent patterns of beach erosion and accretion controlled by ENSO. We nd that approximately one third of all beaches experience signicant erosion during El Niño phases, with the Eastern Pacic particularly vulnerable to widespread erosion (most notably during the large 1997/1998 event). In contrast, La Niña events coincide with signicant accretion for approximately one quarter of all beaches, although conversely drives substantial erosion in south-east Australia and other localized regions. The signicant regional variability in coastal response to ENSO should be considered in light of future projected intensication and shifts in ENSO amplitudes and avors.

map shoreline changes using individual satellite images (typically every 15 days for Landsat), signi cantly increasing the frequency of satellite-derived coastal data.
Here, we use individual satellite images to derive shoreline time-series spanning three Landsat missions  along wave-dominated sandy coastlines around the Paci c Basin. We present a detailed overview of the teleconnections between regional patterns of coastal erosion and accretion in the Paci c and ENSO. To further investigate the mechanisms responsible for the observed regional patterns of shoreline variability, we perform a basin-scale analysis of the variations in coastal wave energy ux and sea-level anomalies during the two opposite phases of ENSO. Lastly, we identify temporal patterns of extreme beach erosion associated with major ENSO events over the past 35 years. To investigate the effect of ENSO on interannual shoreline variability, we evaluated anomalies in seasonally-averaged shoreline position during El Niño and La Niña phases at each transect. The Multivariate ENSO Index (MEI), considered as one of the most complete indices describing ENSO 30 , was used to identify El Niño and La Niña periods during the past 35 years (Methods). We performed a statistical test (non-parametric Wilcoxon signed rank test) to assess whether shoreline anomalies during each phase were signi cantly different (at the 95% con dence level) to the long-term average. Also, with strong seasonal shoreline variability evident particularly in the northern hemisphere, where energetic wave conditions primarily occur during the boreal winter (DJF) 31 , the anomaly analyses were carried out for both boreal winter shorelines as well as for all four seasons (see Methods   1 summarizes the regional sandy coastline response to El Niño and La Niña phases for both the boreal winter-only analysis as well as all four seasons. Considering the entire Paci c, 69% of all transects indicate beach erosion during El Niño phases, of which 33% are statistically signi cant in both analyses (i.e., boreal winter and all seasons). Focusing on the regional variability, the Eastern Paci c (USA, Mexico, Peru and Chile) in particular is identi ed as showing a particularly coherent erosion response to El Niño phases, with 50% of transects (boreal winter) showing signi cant erosion (Fig. 1a). Concurrently, broad accretion is observed during El Niño phases along sandy beaches in south-east Australia, where more than 70% of transects (all seasons) experience signi cant beach accretion (Fig. 1c).
By contrast, during La Niña phases, 60% of all transects indicate beach accretion, of which 26% (boreal winter) and 24% (all seasons) are statistically signi cant. Positive shoreline anomalies are predominantly found in the Eastern Paci c, while widespread erosion is observed during La Niña phases along southeast Australia (45% of transects experiencing signi cant erosion for all seasons, Fig. 1d). Along the coasts of smaller island nations where suitable Landsat imagery were available (New Zealand, Japan and Hawaii), the shoreline data shows a mixed response to ENSO. While the proportion of statistically signi cant anomalies is not su cient to infer a coherent regional response at these islands, it can be observed that Hawaii tends to erode during El Niño and accrete during La Niña phases (all seasons), as does Japan.
Changes in oceanographic forcing associated with ENSO As sandy beaches generally erode in response to energetic waves and/or elevated water levels 7,32,33 and accrete during calmer periods 34,35 , the spatial distribution of these two key oceanographic variables were analysed at the same coastal locations as the satellite shoreline analysis. This was undertaken using ECMWF global products, the ERA5 wave reanalysis and daily sea-level anomalies from satellite altimetry (detrended to reduce sea-level rise, see Methods). Figure 2 shows regional differences in wave energy ux during positive and negative ENSO phases, considering both boreal winter and all seasons. The impact of ENSO on sea-level anomalies across Paci c regions is indicated in Figure 3. The widespread regional erosion during El Niño phases observed particularly in the Eastern Paci c is consistent with a concurrent increase in wave energy ux in USA and Mexico (22% and 14% increase during boreal winter, respectively) and signi cant increases (+5-10 cm) in sea-level anomalies along the whole Eastern boundary. In Peru and Chile, wave energy ux shows a mixed response to El Niño, suggesting elevated sea levels may control the observed erosion response in these regions during El Niño phases. This likely accounts for the lesser widespread shoreline erosion observed in Peru and Chile (47% and 33% of all transects, all seasons) compared to USA and Mexico (54% and 63%, boreal winter). During La Niña phases, deviations in both wave energy ux and sea-level anomalies generally indicate an inverse response in the Eastern Paci c, consistent with the observed overall shoreline accretion during this phase.
In the Western Paci c, sea-level anomalies are relatively weak and for the most part not-statistically signi cant. However, the strong and out-of-phase coastal response to ENSO observed in south-east Australia (Figs. 1c and 1d) is also associated with an increase in wave energy ux (+5%) during La Niña phases and decrease (-5%) during El Niño phases as reported in Figs. 2c and 2d (all seasons). The mixed shoreline response observed along the island nations is likely the result of their widely diverse coastline orientation ( Supplementary Fig. 2) and subsequent exposure to a multitude of wave sources across the Paci c (e.g., New Zealand and Hawaii being exposed to both Northern and Southern Paci c wave climates). This sensitivity to localized wave exposure and associated coastal response is also evident in several sub-regions (e.g., Baja California, northern Peru).

Extreme beach erosion and links to major ENSO events
In addition to regional patterns in coastal erosion and accretion due to ENSO (Fig. 1), shoreline timeseries were analysed to identify extreme beach erosion associated with major ENSO events of the past 35 years. Extreme beach erosion was de ned by the lowest 5% of seasonally-averaged shorelines and collated for each region to highlight widespread extreme erosion in either El Niño, La Niña or Neutral years (Methods). This temporal variability in extreme beach erosion is shown for each region in Figure 4.
Overall, the regional pattern is consistent with the more general shoreline erosion and accretion patterns in Figure 1, with extreme beach erosion associated with major El Niño events in the Eastern Paci c and, inversely, major and prolonged La Niña phases in south-east Australia. In the Eastern Paci c, this analysis reveals that the 1997/1998 event was the most erosive over the past 35 years, during which close to 50% of all transects experienced widespread extreme erosion across the four regions from USA to Chile (Figs. 4e-h). The 2009/2010 and 2015/2016 El Niño events also indicate extreme regional erosion in USA and Mexico, but to a lesser spatial extent (i.e., 29% compared to 42% in USA, and 22% compared to 57% in Mexico for the 2009/2010 and 1997/1998 events, respectively) and are not evident in Peru and Chile.
Also apparent in Figure 4 are examples of extreme coastal erosion in successive years following major ENSO events. For example, in the Eastern Paci c the 1997/1998 El Niño is followed by a La Niña year, but the percentage of extremely eroded transects in 1998/1999 remains high as beaches may still be in a slow recovery period (similarly in Australia, 2012/2013 does not correspond to a major La Niña phase but follows two consecutive years of strong La Niña conditions). Along the complex coastlines of the three smaller island nations, no years in particular stand out in terms of widespread extreme erosion.

Regional variability in ENSO and future Paci c coastal response
We have shown how the El Niño/Southern Oscillation causes large regional variability in coastal erosion and accretion in the Paci c Basin. Interannual shoreline response to ENSO is most coherent along the Eastern Paci c and south-east Australia, but out-of-phase. In view of current ENSO projections, which point towards an increase in the frequency of extreme El Niño events 3

Declarations Data Availability
The full satellite-derived shoreline dataset generated and analysed in the current study is available in the following Zenodo data repository: https://doi.org/10.5281/zenodo.5035998.

Acknowledgements
The authors would like to thank the United States Geological Survey / NASA for providing high-quality open-access data to the scienti c community, Google Earth Engine for facilitating the access to the archive of publicly available satellite imagery, NOAA for maintaining updated time-series of the major climate indices, ECMWF for the reanalysis ERA5 data and multi-mission altimetry dataset, CNES / LEGOS / CLS / AVISO for producing the global tide model FES2014 and Frederic Briol for developing the Python wrapper and the OpenStreetMap project and contributors (https://www.openstreetmap.org) for their extensive geospatial database. The lead author is supported by a UNSW Scientia PhD scholarship. Methods Sandy beaches dataset. A rst pass to identify suitable sandy beaches in the Paci c Basin to investigate regional variability controlled by ENSO was undertaken based on two primary criteria: (1) only wavedominated sandy coastlines (i.e., coastlines primarily controlled by wave, rather than tidal, processes); and (2) Figure 1d. For example, islands in the south-west Paci c (including Papua New Guinea and the Philippines) and in the Artic, as well as the shores of Central America, do not contain su cient images (minimum set at 600) to achieve continuous monitoring over the 35-year study period. Therefore, these regions were also discarded from the analysis. Based on the two criteria described above, a total of 7 regions were identi ed: south-east Australia, New Zealand, Chile, Peru, Mexico's west coast, USA's west coast and Japan. For each region, a dataset of sandy beach locations was compiled using the OpenStreetMap database 47 and manually digitized where necessary. The sandy shoreline locations were then used to generate 100 m alongshore-spaced shore-normal transects.
From a rst-pass dataset of more than 100,000 beach transects using the two criteria above, a secondary, more-detailed, ltering step was applied based on the following conditions: Transects where shoreline time-series (see next section on shoreline mapping) contained insu cient data for beach slope estimation (and hence tidal correction). This corresponded to a de ned lower limit of 180 shorelines at each transect, or 5 shorelines/year [75% of the total ltered transects in this second pass] Transects where the shoreline time-series indicated a strong linear trend of erosion or accretion (threshold ±3m/year based on 1 ), that are likely evidence of non-ENSO in uences on the shoreline (e.g. major beach nourishments, tectonic shifts, land subsidence) [10% of ltered transects] Transects along small pocket beaches (less than 500 m in length) that may be dominated by cellular beach circulation effects [5% of ltered transects] Transects along high-energy meso-macrotidal beaches (but still classi ed as wave-dominated) as these require a more sophisticated water level correction for accurate shoreline mapping 48 . These include sandy beaches along the Oregon and Washington coasts as well as localised stretches of coastline along New Zealand's West coast [10% of ltered transects] This amounted to a nal dataset of 72,000 beach transects across the 7 regions. Supplementary Figure 2 illustrates the spatial extent of this dataset, including the number of transects and their orientation for each region.
Satellite-derived shorelines. The CoastSat open-source Python toolkit 29 (available at https://github.com/kvos/CoastSat) was used to download Landsat images between 1984 and 2020 at each beach (using the Google Earth Engine API 23 ) and automatically map the position of the shoreline on each image. Brie y, CoastSat's detection algorithm combines a supervised image classi cation and a sub-pixel resolution border segmentation to map the position the instantaneous boundary between sand and water with a shoreline accuracy of 10-15 m 26 . A noteworthy improvement over the default CoastSat toolbox was to create a more generalised classi er (originally trained on 5 Australian beaches) using tens of images from sites across the 7 regions, including Chilean black sand beaches, bright sand Atacama Desert beaches, Mexican tropical beaches as well as Japanese Artic beaches. Thus, this extended classi er was capable of identifying sandy pixels at a wide range of coastal environments, making the shoreline detection more robust at the Paci c-basin scale. The intersection of the two-dimensional shorelines with the shore-normal transects were subsequently calculated to produce time-series of crossshore shoreline change at each transect.
Tidal correction. Since every satellite image was taken at a different stage of the tide, the time-series of shoreline change obtained in the previous step were tidally-corrected to remove high-frequency shoreline changes related to tidal variations. To do this, water levels at the time of image acquisition and an estimate of the beach slope are necessary. Since obtaining measured water levels for each beach at the basin scale is not plausible, modelled tides from the FES2014 global tide model 46 were used. The beach slope along each transect was then estimated from the satellite-derived shorelines and corresponding tide levels using a method described in 28 (available at https://github.com/kvos/CoastSat.slope). This technique applies a frequency domain analysis to the shoreline time-series to nd the optimum slope that, when used for tidal correction, minimises high-frequency tidal uctuations relative to lowerfrequency erosion/accretion signals.
Anomalies in shoreline position during ENSO phases. To evaluate the effect of ENSO on shoreline changes, the seasonally-averaged shoreline positions during El Niño and La Niña phases were compared to the long-term average with a distribution-free statistical test (Fig. 1). Two separate anomaly analyses were performed, one including shorelines from all four seasons and one using only the boreal winter (DJF) averages to address the strong seasonality in shoreline variability of the northern hemisphere 31 . El Niño and La Niña phases were de ned by the periods of time during which the intensity of the Multivariate ENSO Index 30 (MEI) was larger than half of its standard deviation (similar to 42 ). Supplementary Figures 3 and 4 show examples of detrended time-series of shoreline change along two different transects located respectively along Ocean Beach (San Francisco, USA) and Narrabeen (Sydney, Australia). The MEI time-series and the observed anomalies in shoreline position are also illustrated in the gures.
These two examples show that at a seasonal northern hemisphere beach like Ocean Beach ( Supplementary Fig. 3), the ENSO anomalies are more pronounced during the boreal winter (DJF), while at a southern hemisphere beach with no apparent seasonality like Narrabeen ( Supplementary Fig. 4), the anomalies are evident during all four seasons.
To test if the observed anomalies are statistically different than the long-term average, the shoreline positions during each ENSO period were compared to the long-term average with a Wilcoxon signed-rank test -the distribution-free equivalent of a one-sample t-test. The resulting p-value indicate the probability of the mean of the population being greater or less than (one-tailed) the long-term average. The statistical test was repeated for both ENSO phases and evaluated for the boreal winter (DJF) anomalies as well as anomalies from all four seasons, generating four p-values for each transect, as presented in the four subplots of Figure 1 in the main manuscript.
Anomalies in wave energy ux during ENSO phases. A similar analysis was performed to quantify the effect of ENSO on wave energy ux around the Paci c (Fig. 2). Six-hourly time-series of signi cant wave height and mean wave period from the 50 km-resolution ECMWF ERA5 reanalysis 45 were used to calculate wave energy ux in watts per metre of wave front: where kg/m 3 , g is the gravitational constant, signi cant wave height and mean wave period. An example of wave energy ux time-series at the closest grid point to Ocean Beach (San Francisco, USA) is shown in Supplementary Figure 5. The statistical signi cance of the observed deviations from the longterm average wave energy ux during both ENSO phases was also quanti ed with a Wilcoxon signedrank test as described above.
Sea-level anomalies during ENSO phases. Time-series of daily sea-level anomalies based on satellite altimetry provided by ECMWF and Copernicus Climate Change Service between 1993 and 2020 were used to investigate corresponding sea level anomalies during ENSO phases. Sea-level anomaly time-series were detrended to reduce any sea-level rise signal and the differences between El Niño and La Niña phases were computed on the subsequent seasonally-averaged time-series. An example of sea-level anomaly time-series at the closest grid point to Ocean Beach (San Francisco, USA) is illustrated in Supplementary Figure 6. As per above, the statistical signi cance of the observed differences in sea level anomalies was also evaluated with a Wilcoxon signed-rank test.
Regional extreme beach erosion. Seasonally-averaged time-series of shoreline variability along each transect were rst divided into yearly bins, de ned from 1 August to 31 July. For each yearly bin, the most eroded datapoint was selected and compared to the lowest 1% and 5% shoreline positions in the record (Fig. 4). To classify El Niño, La Niña or Neutral years, the same approach as in 49 was used and selected the years when the Oct-Nov-Dec (OND) average of the MEI was above (for El Niño years) or below (for La Niña years) half of its standard deviation. Figure 1 Regional patterns of shoreline response to ENSO in the Paci c. Teleconnections between shoreline variability and ENSO along wave-dominated sandy beaches in the Paci c during the boreal winter season only (a,b) and during all four seasons (c,d). The pie charts indicate the proportion of transects that experience enhanced erosion (red) or accretion (blue) during El Niño (a,c) and La Niña (b,d) phases for each region. The proportion of transects along which the anomalies are statistically signi cant at a 95% con dence level (Methods) are hashed.

Figure 2
Teleconnections between wave energy ux and ENSO phases in the Paci c. The maps show the percentage change in mean wave energy ux during the boreal winter season only (a,b) and during all four seasons (c,d) for El Niño (a,c) and La Niña (b,d) phases. Anomalies are calculated along a 50 kmgrid using ERA5 wave data between 1979-2020 (Methods). The anomalies that are statistically signi cant at a 95% con dence interval are circled. A probability distribution function (PDF) of all the grid points belonging to each region is plotted as well, the colour of the distribution represents the position of the mean along the colour bar (e.g., grey if the mean of the distribution is between -5 and 5% anomaly).

Figure 3
Teleconnections between sea-level anomalies and ENSO phases in the Paci c. The maps show differences in detrended sea-level anomalies during the boreal winter season only (a,b) and during all four seasons (c,d) for El Niño (a,c) and La Niña (b,d) phases calculated along a 50 km-grid using the ECMWF multi-mission altimetry dataset between 1993-2020 (Methods). The anomalies that are statistically signi cant at a 95% con dence interval are circled. A probability distribution function (PDF) of all the grid points belonging to each region is plotted as well, the colour of the distribution represents the position of the mean along the colour bar (e.g., grey if the mean is between -1 and 1 cm anomaly).