Pacific shoreline erosion and accretion patterns controlled by El Niño/Southern Oscillation

In the Pacific Basin, the El Niño/Southern Oscillation (ENSO) is the dominant mode of interannual climate variability, driving substantial changes in oceanographic forcing and impacting Pacific 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 38 years of Landsat imagery to map shoreline variability around the Pacific Rim and identify coherent patterns of beach erosion and accretion controlled by ENSO. On the basis of more than 83,000 beach transects covering 8,300 km of sandy coastline, we find that approximately one-third of all transects experience significant erosion during El Niño phases. The Eastern Pacific is particularly vulnerable to widespread erosion, most notably during the large 1997/1998 El Niño event. By contrast, La Niña events coincide with significant accretion for approximately one-quarter of all transects, although substantial erosion is observed in southeast Australia and other localized regions. The observed regional variability in the coastal response to ENSO has important implications for coastal planning and adaptation measures across the Pacific, particularly in light of projected future changes in ENSO amplitude and flavour. The El Niño/Southern Oscillation drives coherent patterns of beach erosion and accretion around the Pacific Rim, according to analysis of satellite imagery covering over 8,300 km of sandy coastline.

In the Pacific Basin, the El Niño/Southern Oscillation (ENSO) is the dominant mode of interannual climate variability, driving substantial changes in oceanographic forcing and impacting Pacific 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 38 years of Landsat imagery to map shoreline variability around the Pacific Rim and identify coherent patterns of beach erosion and accretion controlled by ENSO. On the basis of more than 83,000 beach transects covering 8,300 km of sandy coastline, we find that approximately one-third of all transects experience significant erosion during El Niño phases. The Eastern Pacific is particularly vulnerable to widespread erosion, most notably during the large 1997/1998 El Niño event. By contrast, La Niña events coincide with significant accretion for approximately one-quarter of all transects, although substantial erosion is observed in southeast Australia and other localized regions. The observed regional variability in the coastal response to ENSO has important implications for coastal planning and adaptation measures across the Pacific, particularly in light of projected future changes in ENSO amplitude and flavour.
Sandy coasts are estimated to comprise 31% of coastal environments worldwide 1 , of which the majority are classified wave dominated 2 . These coasts are particularly vulnerable to fluctuations in ocean wave energy and water levels, which drive cycles of erosion and accretion at episodic, seasonal, interannual and decadal timescales, impacting adjacent infrastructure and beach habitats. The interannual timescale is of particular interest as it is closely linked to Earth's climate and internal modes of climate variability. In a changing climate, a likely change in pattern of these important climate drivers 3,4 , coupled with projected changes in storminess 5,6 and rising sea levels, could exacerbate coastal erosion 7 and threaten the future resilience of many coastal communities worldwide 8,9 . In the Pacific Basin, El Niño/Southern Oscillation (ENSO) is the dominant mode of interannual climate variability and has teleconnections with a broad range of atmospheric and oceanic processes along coastal regions 10 , influencing nearshore wave climates 11 , sea-level anomalies 12 and river discharge 13 . Yet our understanding of how sandy coasts in the Pacific respond to these basin-scale changes has to date been limited to a relatively small number (~50 open-coast beaches) of long-term beach monitoring sites, located predominantly along developed coasts in North America, Australia and Japan [14][15][16][17][18][19][20][21][22] , with no observational coverage in Central and South America. Recent innovations in cloud data platforms 23 and remote-sensing algorithms [24][25][26] have opened the possibility for coastal change to be quantified at unprecedented global scales using satellite imagery. While this approach has been successfully used to identify global long-term trends in shoreline change over several decades 1,27 , a major factor that has to date limited the temporal resolution of satellite-derived coastal Article https://doi.org/10.1038/s41561-022-01117-8

Regional patterns of shoreline response to ENSO
Wave-dominated sandy beaches across the Pacific were identified and mapped on the basis of the availability of Landsat images and use of the relative tidal range (RTR = spring tidal range/average wave height) to distinguish wave-dominated (RTR < 3) from more tide-dominated coasts (Methods and Supplementary Fig. 2.1). The CoastSat toolbox 30 was then used to automatically map the shoreline position on Landsat 5 (1984Landsat 5 ( -2013, Landsat 7 (1999Landsat 7 ( -2022 and Landsat 8 (2013Landsat 8 ( -2022 images. A validation of the satellite-derived shorelines against in situ surveys around the Pacific Basin is provided in Supplementary Information. The resulting time series of shoreline change were tidally corrected along 100-m-spaced cross-shore transects using a global tide model and a satellite-derived estimate of the average local beach slope 28 , amounting to 83,677 beach transects in total across the Pacific (equivalent to 8,367 km of sandy coast; full dataset shown in Supplementary Fig. 2.2).
To investigate the effect of ENSO on interannual shoreline variability, we evaluated anomalies in shoreline position during El Niño and La Niña phases at each transect. The Multivariate ENSO Index (MEI), considered one of the most complete indices describing ENSO 31 observations (including interannual variability) has been the effect of high-frequency tidal fluctuations, requiring the use of satellite composite images averaged over large temporal windows (annually). However, new methods for beach slope estimation and tidal correction at the global scale 28 now make it possible to map shoreline changes using individual satellite images (typically every 16 days for each Landsat satellite), significantly increasing the frequency of satellite-derived coastal data 29 .
In this Article, we use individual satellite images to derive shoreline time series spanning three Landsat missions (1984-2022) along wave-dominated sandy coastlines around the Pacific Basin. We present a detailed overview of the teleconnections between regional patterns of coastal erosion and accretion in the Pacific 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 flux and sea-level anomalies during the two opposite phases of ENSO. Last, we identify temporal patterns of extreme beach erosion associated with major ENSO events over the past 38 years.  32 , anomaly analyses were carried out for boreal winter shorelines as well as shorelines for all four seasons (Methods and Supplementary Information 3.2). Figure 1 summarizes the regional sandy coastline response to El Niño and La Niña phases for the boreal winter-only shorelines as well as shorelines from all four seasons. Considering the entire Pacific, 70% of all transects indicate beach erosion during El Niño phases, of which 33% are statistically significant in the two analyses (boreal winter and all seasons). Focusing on the regional variability, the Eastern Pacific (United States, Mexico, Peru and Chile) is identified as showing a particularly coherent erosion response to El Niño phases, with 50% of transects (boreal winter) showing significant erosion (Fig. 1a). Concurrently, broad accretion is observed during El Niño phases along sandy beaches in southeast Australia, where 75% of transects (all seasons) experience significant beach accretion (Fig. 1c).
By contrast, during La Niña phases, 60% of all transects indicate beach accretion, of which 25% are statistically significant. Positive shoreline anomalies are found predominantly in the Eastern Pacific, while widespread erosion is observed during La Niña phases along southeast Australia (48% of transects experiencing significant erosion for all seasons; Fig. 1d). Along the coasts of smaller island nations where suitable Landsat imagery was available (New Zealand, Japan and   are statistically significant at a 5% significance level are outlined with a circle. A probability distribution function 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 (for example, grey if the mean of the distribution is between -5% and 5% anomaly).
Article https://doi.org/10.1038/s41561-022-01117-8 Hawaii), the shoreline data show a mixed response to ENSO. While the proportion of statistically significant anomalies is not sufficient to infer a coherent regional response at these islands, it can be observed that Hawaii tends to erode during El Niño phases 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,33,34 and accrete during calmer periods 35,36 , the spatial distributions of these two key oceanographic variables were analysed at the same coastal locations as the satellite shoreline analysis. This was undertaken using the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth-generation wave reanalysis (ERA5) and daily sea-level anomalies from satellite altimetry observations (detrended to reduce sea-level rise; Methods), with both datasets provided at a global scale. Figure 2 shows regional differences in wave energy flux during positive and negative ENSO phases, considering both boreal winter and all seasons. The impact of ENSO on sea-level anomalies across Pacific regions is indicated in Fig. 3. The widespread regional erosion during El Niño phases observed particularly in the Eastern Pacific is consistent with a concurrent increase in wave energy flux in the United States and Mexico (22% and 14% increases during boreal winter, respectively) and significant increases (+5 to +10 cm) in sea-level anomalies along the whole eastern boundary. In Peru and Chile, wave energy flux shows a weak response to El Niño. This probably accounts for the lesser widespread shoreline erosion observed in Peru and Chile (44% and 33% of all transects, all seasons) compared with the United States and Mexico (54% and 67%, boreal winter). During La Niña phases, deviations in both wave energy flux and sea-level anomalies generally indicate an inverse response   (Methods). The anomalies that are statistically significant at a 5% significance level are outlined with a circle. A probability distribution function 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 (for example, grey if the mean is between -1 and 1 cm anomaly).
In the Western Pacific, sea-level anomalies are for the most part not statistically significant. However, the strong and out-of-phase shoreline response to ENSO observed in southeast Australia (Fig. 1c,d) is also associated with an increase in wave energy flux (+7%) during La Niña phases and a decrease (-6%) during El Niño phases as reported in Fig.  2c,d (all seasons). The mixed shoreline response observed along the island nations is probably the result of their widely diverse coastline orientation ( Supplementary Fig. 2.2) and subsequent exposure to a multitude of wave sources across the Pacific (for example, New Zealand and Hawaii being exposed to both northern and southern Pacific wave climates). This sensitivity to localized wave exposure and associated coastal response are also evident in several sub-regions (for example, Baja California and 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 time series were analysed to identify extreme beach erosion associated with major ENSO events of the past 38 years. Extreme beach erosion was defined by the lowest 5% of seasonally   Year Year 1984 1987 1990 1993 1996 1999 2002 2005 2008 2011 2014 2017 2020 1984 1987 1990 1993 1996 1999  Article https://doi.org/10.1038/s41561-022-01117-8 averaged shorelines and collated for each region to highlight widespread extreme erosion during El Niño, La Niña or neutral years (Methods). This temporal variability in extreme beach erosion is shown for each region in Fig. 4. Overall, the regional pattern is consistent with the more general shoreline erosion and accretion patterns in Fig. 1, with extreme beach erosion associated with major El Niño events in the Eastern Pacific and, inversely, major and prolonged La Niña phases in southeast Australia. In the Eastern Pacific, this analysis reveals that the 1997/1998 Eastern Pacific El Niño event was the most erosive over the past 38 years, during which close to 50% of all transects experienced widespread extreme erosion across the four regions from the United States to Chile (Fig.  4e-h). The 2009/2010 Central Pacific and 2015/2016 mixed (Eastern Pacific/Central Pacific) El Niño events also indicate widespread regional erosion in the United States and Mexico, but to a lesser spatial extent than the 1997/1998 Eastern Pacific event. On the basis of these observations, in the United States and Mexico, 43% and 57% of the transects, respectively, experienced extreme erosion (defined as the lowest seasonally averaged shoreline being below the 5th percentile) during the 1997/1998 event. By comparison, during the 2009/2010 and 2015/2016 events, this number amounts to only 29% and 18%, respectively, for the United States and 22% and 23% for Mexico. In Peru and Chile, the spatial extent of the extreme erosion observed during the 1997/1998 event (~50% of transects) has not been observed since then, with the 2009/2010 and 2015/2016 events not resulting in particularly widespread erosion.
Also apparent in Fig. 4 is the role of antecedent conditions on extreme erosion, with observations of persistent erosion in years following major ENSO events. For example, in the Eastern Pacific, the 1997/1998 El Niño is followed by a La Niña year, but the percentage of extremely eroded transects in that year remains high as beaches may still be in a slow recovery period. Similarly in Australia, the most widespread erosion (47% of transects experiencing extreme erosion) is observed in 2012/2013, which does not correspond to a major La Niña phase but follows two consecutive years of strong La Niña conditions. This highlights how El Niño and La Niña can trigger prolonged erosion phases on sandy coastlines, as has been observed at local sites using long-term in situ data 37 . Future work is needed to clarify the complex lagged responses among oceanic forcing, erosion and longer-term beach recovery at the Pacific Basin scale.

Regional variability in ENSO and future Pacific coastal response
We have shown how ENSO causes large regional variability in coastal erosion and accretion along the Pacific Rim. Interannual shoreline response to ENSO is most coherent along the Eastern Pacific and southeast Australia, but with opposite phase. In view of current ENSO projections, which point towards an increase in the frequency of extreme El Niño events 3,38 , a shift in their 'flavour' 39 from Eastern Pacific to Central Pacific 40,41 and an increase in frequency of La Niña events 42 , our analyses suggest that the Eastern Pacific and southeast Australia emerge as the sectors of the Pacific Rim most susceptible to enhanced ENSO-driven interannual shoreline variability in a warming climate.
Furthermore, the temporal variability in extreme coastal erosion spanning the past 38 years suggests that, along the Eastern Pacific (California to Chile), recent El Niño events with a strong Central Pacific component 43-46 (the 2003/2004, 2009/2010 and 2015/2016 events) have resulted in less widespread erosion compared with the canonical 1997/1998 Eastern Pacific event. This is in line with the severity of oceanographic forcing observed during each event, also reported in ref. 22 , which shows that during the three strongest events in the study period, the largest increase in wave energy flux and sea-level anomalies along the Eastern Pacific (California to Chile) occurred during the 1997/1998 event, with relatively weaker oceanographic forcings during the 2015/2016 mixed event and substantially weaker forcings during the 2009/2010 Central Pacific event ( Supplementary Fig. 5). Future research should focus on understanding the mechanisms by which different types of El Niño influence oceanographic forcing across the oceanic basin and how these may affect open-coast sandy coastlines in a Central Pacific El Niño-dominated climate.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-022-01117-8.

Sandy beaches dataset
A first pass to identify suitable sandy beaches in the Pacific Basin to investigate regional variability controlled by ENSO was undertaken on the basis of two primary criteria: (1) only wave-dominated sandy coastlines (coastlines controlled primarily by wave, rather than tidal, processes) and (2) the availability of Landsat images for satellite-derived shoreline mapping, due to the inconsistent geographical coverage of the Landsat archive 47 .
Wave-dominated coastlines were classified on the basis of the RTR 48 , calculated here as the ratio between the average deep-water significant wave height, obtained from ERA5 reanalysis 49 and the spring tidal range, obtained from the FES2014 global tide model 50 . Supplementary Fig. 2.1 depicts the spatial distribution of significant wave height ( Supplementary Fig. 2.1a), spring tidal range ( Supplementary  Fig. 2.1b) and the resulting RTR ( Supplementary Fig. 2.1c). Coastlines with an RTR ≤ 3 were classified as wave dominated, while the other coastlines were labelled as tide modified (3 < RTR ≤ 10) or tide dominated (10 < RTR ≤ 50). Regions excluded from this study on the basis of this criterion included northeast Australia, China, Alaska, Canada and Central America.
Regarding image availability, Landsat images suitable for time-series analysis are categorized as Tier 1 due to their consistent quality through time and across instruments. The availability of Tier 1 Landsat images between 1984 and 2020 (Landsat 5, 7, 8) varies across the Pacific Basin, as shown in Supplementary Fig. 2.1d. For example, islands in the southwest Pacific (including Papua New Guinea and the Philippines) and in the Arctic, as well as the shores of Central America, do not contain sufficient images (minimum set at 600) to achieve continuous monitoring over the 38-year study period. Therefore, these regions were also discarded from the analysis. On the basis of the two criteria described in the preceding, a total of seven regions were identified as suitable for this analysis: Australia's southeast coast, New Zealand, Chile, Peru, Mexico's west coast, the west coast of the United States and Japan. For each region, a dataset of sandy beach locations was compiled using the OpenStreetMap database 51 (using the query natural type = 'beach' and surface type = 'sand') and manually digitized where necessary. The sandy shoreline locations were then used to generate 100 m alongshore-spaced shorenormal transects.
From a first-pass dataset of more than 110,000 beach transects using the preceding two criteria, a secondary, more-detailed, filtering step was applied to discard transects on the basis of the following conditions: • Transects where shoreline time series (see next section 'Satellite-derived shorelines') contained insufficient data for beach slope estimation (and hence tidal correction); this corresponded to a defined lower limit of 180 shorelines at each transect, or 5 shorelines per year (75% of the total discarded transects in this second pass) • Transects where the shoreline time series indicated a strong linear trend of erosion or accretion (threshold ±3 m yr -1 on the basis of ref. 1 ), that are likely evidence of non-ENSO influences on the shoreline (for example, major beach nourishments, tectonic shifts and land subsidence) (10% of discarded transects) • Transects along small pocket beaches (less than 500 m long) that may be dominated by cellular beach circulation effects (5% of discarded transects) • Transects along high-energy meso-macrotidal beaches (but still classified as wave dominated) as these require a more sophisticated water-level correction for accurate shoreline mapping 52 ; these include sandy beaches along the Oregon and Washington coasts as well as localized stretches of coastline along New Zealand's west coast (10% of discarded transects) This amounted to a final dataset of 83,677 beach transects (74% of the first-pass dataset) across the seven regions used for the ENSO anomaly analysis. Supplementary Fig. 2.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 30 (available in ref. 53 ) was used to extract Landsat images between 1984 and 2022 at each beach (using the Google Earth Engine Application Programming Interface 23 ) and automatically map the position of the shoreline on each image. Briefly, CoastSat's detection algorithm combines a supervised image classification and a sub-pixel resolution border segmentation to map the position of the instantaneous boundary between sand and water with a shoreline accuracy of 10-15 m (ref. 26 ) (see comparison with in situ shoreline data in the section 'Validation of satellite-derived shorelines'). A noteworthy improvement in this study over the default CoastSat toolbox was to create a more generalized classifier (originally trained on five Australian beaches) using tens of images from sites across the seven regions, including Chilean black-sand beaches, Chilean bright-sand Atacama Desert beaches, Mexican tropical beaches and Japanese snow-covered beaches. Thus, this extended classifier was capable of identifying sandy pixels at a wide range of coastal environments, making the shoreline detection more robust at the Pacific Basin scale. The intersections of the two-dimensional shorelines with the shore-normal transects were subsequently calculated to produce time series of cross-shore 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 50 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 ref. 28 (available in ref. 53 ). This technique applies a frequency domain analysis to the shoreline time series to find the optimum slope that, when used for tidal correction, minimizes high-frequency tidal fluctuations relative to lower-frequency erosion/ accretion signals.

Validation of satellite-derived shorelines
The accuracy of the satellite-derived shorelines obtained with CoastSat has been evaluated in two previous studies 26,52 , and to date, the toolbox has been used in a number of recent publications by other authors to extract shoreline observations along open-coast sandy beaches worldwide [54][55][56][57][58][59][60] . In Supplementary Information 1, we have compiled all the in situ long-term shoreline datasets available along the Pacific Rim: • 16 years at Torrey Pines (San Diego, USA), publicly available and described in ref. 61 • 18 years at Ocean Beach (San Francisco, USA), described in refs. 62,63 • 34 years at Narrabeen-Collaroy (Australia), publicly available and described in ref. 64 • 10 years at Moruya (Australia), described in ref. 65 • 20 years at Tairua Beach (New Zealand), described in refs. 66,67 The comparison between in situ surveyed data and satellite-derived shorelines is presented in Supplementary Fig. 1.1 and indicates a standard deviation error of the order of 10 m across more than 10,000 in situ survey points, consistent with previous validation outside the Pacific Rim 26 . This demonstrates how monthly averaged satellite-derived Nature Geoscience Article https://doi.org/10.1038/s41561-022-01117-8 shorelines are capable of capturing interannual shoreline variability at these five sites. In addition, the range of coastal environments (characterized by tidal range, incident wave energy and beach-face slope) captured by the five validation sites is compared with the variability of the entire analysed dataset in Supplementary Fig. 2.3. This indicates a good coverage of the parameter space, noting that the dataset was restricted to microtidal wave-dominated beaches (see previous section 'Sandy beaches dataset'). It is acknowledged, however, that the validation sites do not include beaches that are optically different (for example, black-sand beaches) or composed of coarse sediment (for example, gravel).

Anomalies in shoreline position during ENSO phases
To evaluate the effect of ENSO on shoreline variability, monthly averaged shoreline positions during El Niño and La Niña phases were compared with the long-term average shoreline position along each transect. El Niño and La Niña phases were defined as the periods during which the intensity of the MEI 31 was stronger than half of its standard deviation, in a similar way to ref. 45 . Monthly averaged shoreline positions were then computed to homogenize the temporal resolution of the satellite-derived time series and avoid a potential bias caused by the high density of shorelines over the more recent part of the record. To test whether the shorelines during the El Niño/La Niña phases are significantly wider or narrower than the long-term average, a statistical test was used, a Wilcoxon signed-rank test-the distribution-free equivalent of a one-sample t test 68 . The resulting P value indicates the probability of the mean of the population being greater or less than (one-tailed) the long-term average. Distributions with a mean greater (less) than the long-term average were labelled as 'erosion' ('accretion'). The statistical test was repeated for both ENSO phases and evaluated for the boreal winter (DJF) anomalies as well as anomalies from all seasons, generating four P values for each transect, as presented in Fig. 1a-d. The two separate anomaly analyses, over DJF-only shorelines and all seasons, were performed to address the strong seasonality in shoreline variability of the Northern Hemisphere 32 and avoid potential biases that may emerge when ENSO is locked in the seasonal cycle (further discussion on this in Supplementary Information 3.2). The statistical robustness of this analysis is discussed in detail and validated in Supplementary Information 4 , while at a Southern Hemisphere beach with no apparent seasonality such as Narrabeen, the anomalies are evident when shorelines from all four seasons are considered. As in situ data are available at these two sites, the same composite analysis is also performed on the in situ record in Supplementary Figs. 3.3.3 and 3.3.4. This analysis demonstrates that the same anomalies and statistical significance found in the satellite-derived time series are obtained using the in situ record, providing further support as to the robustness of the Pacific-wide shoreline analyses.

Anomalies in wave energy flux during ENSO phases
A similar analysis was performed to quantify the effect of ENSO on wave energy flux around the Pacific (Fig. 2). Six-hourly time series of significant wave height and mean wave period (1984-2022) from the 50-km-resolution ERA5 reanalysis 49 were used to calculate wave energy flux in watts per metre of wave front: where ρ = 1,025 kg m -3 , g is the gravitational constant, H s is significant wave height and T is mean wave period. An example of wave energy flux time series at the closest grid point to Ocean Beach is shown in Supplementary Fig. 3.4. To maintain the same timescale as the satellite observations, the six-hourly wave time series were monthly averaged, and the statistical significance of the observed deviations from the long-term average wave energy flux during both ENSO phases was evaluated with a Wilcoxon signed-rank test, as previously described and illustrated step by step in Supplementary Fig. 3.4.

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 2022, extracted on the same grid as the wave data, were used to investigate the influence of ENSO on sea-level anomalies. The time series of sea-level anomalies were detrended to reduce any sea-level-rise signal, and as in the previous analysis, the time series were monthly averaged. An example of sea-level anomaly time series at the closest grid point to Ocean Beach and its step-by-step analysis is provided in Supplementary Fig. 3.5. As previously, the statistical significance 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 first divided into yearly bins, defined from 1 August to 31 July. For each yearly bin, the most eroded datapoint was selected and compared with the lowest 1% and 5% shoreline positions in the record (first and fifth percentiles). The percentage of transects in each of the seven regions that falls below these percentiles is reported in Fig. 4.
To classify El Niño, La Niña and neutral years, the same approach as in ref. 69 was used. The average MEI value during OND (the months over which ENSO is known to peak 70 ) is first calculated and compared with half a standard deviation of the MEI to classify each year into an El Niño year (MEI OND > 0.5σ MEI ), La Niña year (MEI OND < −0.5σ MEI ) or neutral year (|MEI OND | < 0.5σ MEI ). The MEI OND averaged are plotted in Fig. 4i, while El Niño, La Niña and neutral years are colour coded in Fig. 4a-h.

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.4760144. The data are also displayed on an interactive web portal at http://coastsat.wrl.unsw.edu.au/. Source data are provided with this paper.