A Multiscale Approach to Shoreline Prediction

Shorelines respond to a number of “drivers” operating on a variety of time‐scales. For some time‐scales (e.g., seasonal), the driver‐shoreline relationship is often evident; however, at longer time‐scales (e.g., multiannual), the shoreline changes may be superimposed on changes at shorter time‐scales and thus are diﬃcult to identify. Here, we predict shoreline evolution from storm events to decadal time‐scales, using a novel approach based on the Complete Ensemble Empirical Mode Decomposition. This approach identifies and links the primary time‐scales in the model drivers (large‐scale sea level pressure [SLP] and/or waves) with the same time‐scales in the shoreline position. The multiscale approach reproduced shoreline changes at two beaches more skillfully than a common shoreline model when SLP and wave information were used in combination. In addition, the analysis can be applied to climate indices, providing the opportunity to link longer time‐scales with climate patterns (e.g., El Niño Southern Oscillation).

based on SLP to explain wave variability along the coast of Western Europe and extreme beach erosion events. So far, only a few approaches have taken advantage of the link between atmospheric circulation patterns and shoreline changes. For example, Antolınez et al. (2016) proposed a multiscale climate emulator that considers the global atmospheric circulation fields as part of the input of a "one-line" shoreline model. Robinet et al. (2016) developed a statistical model to predict interannual shoreline variability at Truc Vert beach, France, using weather regimes based on the regional atmospheric circulation patterns and showed similar skill to a typical equilibrium model (Yates et al., 2009). Anderson et al. (2018) developed a climate index to relate beach rotation to climate patterns at multidecadal time-scales and showed that large shoreline excursions are related to extreme El Niño winters. More recently, Wiggins et al. (2020) used a combination of the NAO and WEPA index to predict beach rotation (i.e., where one extremity of the beach progradates while the other erodes).
In this study, a new methodology is introduced to predict shoreline changes, in which each time-scale in both drivers and shoreline response is isolated and linked using Complete Ensemble Empirical Mode Decomposition (CEEMD) (Huang et al., 1998;Torres et al., 2011). We used two different model drivers: largescale SLP fields and their gradients, which might be considered a proxy for large-scale climate patterns (Camus et al., 2014a(Camus et al., , 2014b, and bulk wave parameters (H s ,T p ,  ). Our hypothesis is that shoreline position at an individual time-scale can be predicted using the oscillations in the drivers at the same time-scale. The overall shoreline change can be predicted by superimposing changes predicted for each of the relevant time-scales.

Study Sites
Two beaches were considered in this study: Narrabeen Beach, Australia (11 years of shoreline data, Phillips et al., 2017), and Tairua, New Zealand (18 years of shoreline data, Montaño et al., 2020a). Narrabeen Beach is a 3.6 km long embayment located in the Sydney metropolitan area. The beach is micro-tidal (mean spring tidal range = 1.3 m) with coarse sand (D 50 ∼ 0.4 mm). The beach morphology at the profile used in this study (PF6, Turner et al. [2016]) experiences a range of intermediate beach states depending on antecedent wave conditions (Wright & Short, 1984). Shoreline data at Narrabeen were obtained at approximately daily intervals between 2004 and 2015 using an Argus camera system (Harley et al., 2011;Phillips et al., 2017). Tairua Beach is located on the east coast of the North Island of New Zealand. Tairua is a 1.2 km long pocket beach, with D 50 ∼ 0.3 mm that exhibits intermediate beach states (Blossier et al., 2016;Gallop et al., 2011). The beach is micro-tidal with a tidal range varying between 1.2 and 2 m. Shoreline data were obtained using a camera system and the alongshore-averaged cross-shore position was averaged weekly (Blossier et al., 2016(Blossier et al., , 2017Montaño et al., 2020a).

Model Drivers
We considered two model drivers: large-scale two-dimensional SLP fields and gradients from CFSR reanalysis, and wave time-series (H s , T p , and  ). The SLP influence area was identified using the ESTELA method (Pérez et al., 2014). This method evaluates the source and travel time of waves reaching a given location based on the geographic criteria and the two-dimensional wave spectra (IFREMER wave hindcast, Rascle & Ardhuin, 2013). In order to also account for the travel time of swell waves in the analysis, the SLP information is modified according to the isochrones of the average travel time as in Hegermiller et al. (2017) (Figures 1a and 1e). A principal component analysis (PCA) is subsequently performed over the SLP fields and gradients to obtain the principal components-PCs (temporal coefficients; Figures 1c and 1g), associated with the spatial variability patterns, empirical orthogonal functions (EOFs) (Figures 1b and 1f). Here, we use only the first 10 PCs, which explain up to 55% (Narrabeen) and 42% (Tairua) of the overall shoreline variance at these two locations. To fully explain variability, many more PCs are needed, for instance, Rueda et al. (2019) used 568 PCs to explain 95% of the variability in New Zealand. Wave parameters at Narrabeen were obtained from a buoy located 11 km offshore (80 m water depth). The deepwater observations were transformed to 10 m depth contour (immediately offshore of profile PF6) using a SWAN nearshore wave model . Similarly, waves at Tairua were obtained at 10 m water depth using SWAN forced with deepwater waves from WaveWatch III hindcast. For the analysis, all time-series were averaged at a daily scale.

A Model for Shoreline Prediction at Different Time-Scales (SPADS)
SPADS is based on the CEEMD, which is a noise-assisted data analysis based on the Empirical Mode Decomposition (EMD) introduced by Huang et al. (1998). EMD-based methods were designed to identify nonlinear and nonstationary oscillations in data, even with small amplitudes, assuming that simple oscillatory modes of significantly different superimposed frequencies coexist (Huang et al., 1998). EMD decomposes time series into a finite set of "intrinsic mode functions" (IMFs), representing different time-scales with varying amplitudes and frequencies. These are similar to spectral decomposition, but without the requirement of stationarity. The last IMF is a monotonic function, where it is no longer possible to extract further features from the time-series, and is generally treated as a trend. This means that the method does not require detrending of the time-series. EMDs have shown good performance when compared with established methods used in geophysical and climatology studies (Carmona & Poveda, 2014;Ji et al., 2014;Wu et al., 2008). CEEMD improves on EMD, by assist-MONTAÑO ET AL.  ing the IMF separation through the addition of Gaussian white noise to the signal and the true IMF is calculated as the mean of an ensemble of tests (100 in our study) (Torres et al., 2011;Wu & Huang, 2009).
We used white noise amplitudes between 0.1 and 0.5 of the signal standard deviation. Then, a test that identifies the statistically significant oscillations (Wu & Huang, 2004) was performed, and only the IMFs that satisfied a significance level higher than 95% were selected for further analysis. To summarize, the methodology consists of five steps (see Figure S1): 1. Decompose the time-series of the shoreline ( Figure 2) and drivers-SLP and/or wave informationusing CEEMD, in order to isolate time-series representing the different time-scales. Then, apply a statistical test (Wu & Huang, 2004) to identify the significant IMFs (see Step 1 in flowchart, Figure S1).

MONTAÑO ET AL.
10.1029/2020GL090587 4 of 11 2. Separately reconstruct the shoreline oscillation (S IMF,j ) at each significant time-scale, using the drivers (Y IMF,i ) with similar time-scale ( Figure 3). The normalized IMFs of the drivers are used to find the coefficients c that best fit the shoreline position at a specific time-scale, that is, IMF,   c) and (e)-(g) Shoreline oscillation (pink line, left axis) and an example of a driver with a strong influence (blue line, right axis). Dotted black line shows the averaged shoreline prediction using the SLP and wave information. Pink shades represent the standard deviation of the white noise runs. T s is the period (days) of the shoreline oscillation. (d) and (h) behavior of the PC1 (black line, right axis) and SOI (color bars, left axis). Cross-correlation between the driver and the shoreline (R SD ) and correlation coefficient (R) are displayed within each panel.  IMF related to the residual (pink line) explains up to 46% of the shoreline signal. Seasonal (∼181 days), annual (∼363 days), and bi-annual (∼709 days) oscillations are observed in the shoreline data, and these oscillations explain about 11%, 18%, and 12% of the shoreline variability, respectively (Figures 2b-2d). The smaller %E.V of these individual time-scales are due to the residual IMF explaining the larger component of the overall shoreline variance.

Shoreline Position at Individual Time-Scales
We hypothesize that shoreline changes at a specific time-scale can be predicted using the drivers with a similar temporal oscillation. For instance, seasonal shoreline changes might be predicted using only seasonal oscillations in the drivers while the total shoreline change results from the summation of the different time-scales (e.g., seasonal, annual, bi-annual, and decadal). Figures 3a-3c and 3e-3g show the individual shoreline oscillations (pink line, IMF-S) for a specific time-scale (same as Figure 2) and the reconstructed shoreline (dotted black line), using both the time-series of the SLP fields and gradients (i.e., PCs Figures 1c and 1g) and wave bulk parameters (Figures 1d and 1h). Some oscillations, for instance, 980 days at Tairua (Figure 3h), are more difficult to reconstruct since only a few drivers are significant at that scale and, opposite to the annual oscillations, their signal shows stronger sensitivity to the white noise selection (pink shade). Blue lines in Figure 3 (right axis) show some of the drivers that had a strong cross-correlation coefficient (R SD ) over the time with the shoreline position at the individual time-scale, which is reflected in the coefficients (c i , Step 3) when reconstructing the shoreline. Moreover, changes in the amplitude of the drivers (blue lines) are reflected in the shoreline position (pink lines), allowing us to account for the nonstationarity of the wave climate (e.g., Figure 3c). In Figures 3a-3c, 3e-3g, the driver (blue line) and subsequent shoreline response (pink line) are not in phase, and the largest cross-correlations (R SD ) are found at time-lags different than zero days. In Section 3.5, the relationship of some of the drivers with the Southern Oscillation Index (SOI) index at long time-scales is discussed (Figures 3d and  3h).

Hindcasting the Total Shoreline Change
SPADS reproduces the overall shoreline position as the sum of the shoreline reconstructed at each timescale. Figures 4a and 4f show the predicted shoreline changes at Narrabeen and Tairua using SLP and wave information as the basis for prediction. We repeated the analysis using: 1) SLP information only (Figures 4b  and 4g); and 2) wave parameters only (H s , T p , θ) (Figures 4c and 4h). The new approach is compared with the shoreline model ShoreFor , which is applied at coastal environments globally (Dodet et al., 2019;Montaño et al., 2020a;Splinter et al., 2014) due to its simplicity (it requires only H s , T p , sediment size, and only uses three free coefficients) and excellent performance at reproducing shoreline change. At Narrabeen, our model performed better than ShoreFor when SLP information was included (Figures 4a and  4b). When only waves were used at Narrabeen, the model performance was poorer. The model was able to reproduce some accretion (e.g., 2006) and erosion events (e.g., 2007) when SLP information was used. The accretion/erosion events were underestimated by ShoreFor and by our model using only waves as a driver.
For Tairua the best model result was also obtained when SLP and waves were used (Figure 4f), but, contrary to Narrabeen, similar results were obtained when the model used only SLP (Figure 4g) or only wave information (Figure 4h). At Tairua, ShoreFor performed similar to the only-SLP and only-wave cases. Quantile analysis has been performed as it highlights performance for extreme events (Figures 4e and 4j, left panels), showing that our approach captured these events better (a perfect model follows the dashed black line).

Predictive Capability
We also tested the model prediction capability (gray area, Figures 4d and 4i) using SLP and waves. The calibration data set of the models was reduced and the last 2.5 years of each time-series were used to test the model predictions. The last years of shoreline change at Narrabeen show events not present in the calibration period (large accretion-erosion sequence in 2014-2015), making predictions particularly challenging. Based on common metrics (e.g., R 2 , RMSE), at Narrabeen our model had a better performance during the forecast period than ShoreFor, but the contrary was observed at Tairua. Nonetheless, our model captures the extremes more accurately (closer to the black dashed line, Figures 4e and 4j, right panels). We suggest that one of the reasons for this is that ShoreFor tries to optimize a parameter for all time-scales, whereas this approach treats each time-scale separately and is hence able to better capture both short and long time-scales. In addition, SPADS could better reproduce shoreline change (including the extremes) when SLP information is considered (Figures 4a and 4b), since, some processes like surge or more complex wave information can be captured by SLP and help to more accurately predict erosion.
Changes in seasonal and interannual wave variability are expected under future wave climate scenarios (Morim et al., 2019), which are likely to affect the predictive capability of shoreline models (D'Anna et al., 2020;Splinter et al., 2016). Since our model identifies several time-scales, changes in the dominant time-scale can still be captured, but only if the time-scale was present during the calibration period. Additionally, our approach does not depend on the preexisting shoreline position, preventing error propagation. Nonetheless, similar to equilibrium models , our model is sensitive to the length of the data used in the calibration particularly for the longer time-scales.

Extracting Information at Different Time-Scales
Our results show that model performance improved when SLP and wave information are both used. The two model drivers contain related information but there are also differences, which might explain why their combined use improves model skill. The SLP might contain information of fluctuations related to the mean sea level (Robinet et al., 2016;Ruggiero et al., 2001;Serafin & Ruggiero, 2014), not present in the wave data. Even though both study sites do not experience particularly large storm surge, sea level variations can affect erosional patterns (Coco et al., 2014). We expect this information to be even more important at other locations. Moreover, the PCs-EOFs obtained from the SLP track wave information from large-scale generation areas and contain information of large-scale climate patterns (e.g., Camus et al., 2014b;Pérez et al., 2014;Rueda et al., 2019) Also, nearshore wave information (e.g., 10 m water depth) reflects local bathymetric effects that are not captured by the large-scale SLP. Conversely, bulk wave parameters oversimplify other factors like, for example, the possible bi-modality of the wave spectra (e.g., Montaño et al., 2020b;Wiggins et al., 2020). Overall, both drivers are complementary.
We analyzed the influence of the two model drivers in the shoreline prediction at specific time-scales. For instance, at Narrabeen, the largest cross-correlations (R SD ) between drivers and changes in the shoreline position are found at time-scales shorter than annual: ∼45, 85, and 181 days for H s (R SD = −0.47; −0.41; −0.44, Figure 3a); while very low R SD were found with the PCs (SLP fields and gradients) at these time-scales. Narrabeen has been identified as a storm dominated beach (Splinter et al., 2014;Turner et al., 2016), with a weak annual signal in the wave height ( Figure 1d). Our analysis shows that the annual scale is better captured by atmospheric patterns (PC1, Figure 3b) with a high R SD = −0.75 compared to H s (R SD = −0.31). At the bi-annual time-scale we found a significant R SD with PC1 (−0.37) while H s showed a similar nonstationary pattern than the shoreline with R SD = −0.76 (Figure 3c). This interannual erosion/accretion patterns (∼2 years and longer) have been attributed to a more energetic wave climate during la Niña years (Ranasinghe et al., 2004;Turner et al., 2016).
Similar to Narrabeen, Tairua also displayed largest R SD with H s at the short time-scales, without a significant contribution from SLP, showing that SLP information is more important at longer time-scales. At Tairua, contrary to Narrabeen, both drivers displayed a good R SD = 0.60 (PC1), −0.63 (H s , Figure 3f) at the annual scale, which is not unexpected since Tairua displays a clear annual cycle in both waves and shoreline changes. Although the short-term oscillations improve the predictions and better capture some extremes, most of the E.V of the shoreline at both beaches are associated to longer time-scales. Thus, the weak correlation between H s and the shoreline at the annual scale in Narrabeen is reflected in a poor prediction when only wave information is used, while at Tairua the overall shoreline response can be well predicted using only SLP or only wave information.
Previous models Splinter et al., 2014) also extract information on the dominant time-scale of shoreline change using the so-called memory decay parameter (φ), which reflects the shoreline response to wave conditions in terms of sediment exchanges with the surf zone. At Narrabeen φ is about 15-30 days, which can be related to the frequency of the storms (Davidson et al, , 2017Phillips et al., 2017;Splinter et al., 2014), while at Tairua the φ was 1 order the magnitude higher (220 days). It is striking that SPADS, which does not explicitly use a beach memory parameter, captures a relationship between wave parameters and shoreline response at similar time-scales.
Different from previous studies, we analyzed the influence of the ENSO (through the SOI) at individual time-scales. For instance, the 709 days oscillation in the shoreline was significantly correlated with the PC1-EOF1 (R SD = −0.37), which in turn had a significant correlation coefficient (R = 0.46, p < 0.001) with the SOI at this time-scale ( Figure 3d). Although some of the PCs-EOFs have a small E.V (<5%), they can still improve shoreline predictions at some specific time-scale. For instance, the 3.5 years oscillation (1,291 days) had a strong R SD with the PCs-EOFs 5, 6 and 8 (0.69, 0.59 and −0.46) which in turn had a significant R with the SOI (−0.48, −0.4, and 0.55) (supporting information).
In New Zealand, the influence of the ENSO has been shown to be relevant (Godoi et al., 2016;Rueda et al., 2019). We found that the SOI has a significant correlation at the time-scale of ∼980 days with the PC-EOF1 (R = 0.47) (Figures 3c and 3d). In addition, the R SD at this time-scale was large for T p (0.77) which was also significantly correlated with the SOI (supporting information).

Conclusions
We propose a robust modeling approach to shoreline prediction (SPADS) that bridges previous modeling efforts to predict short-term shoreline changes driven by waves and long-term changes driven by largescale atmospheric patterns. SPADS identifies characteristic oscillations in the shoreline and associated drivers, predicts shoreline changes at individual time-scales, and allows us to identify the drivers with the largest contribution. SPADS can make use of different drivers, with the best model performance obtained when SLP and bulk wave information are used. The focus on resolving and predicting individual time-scales facilitates the exploration of oscillations longer than annual, which have a remarkable effect on the overall shoreline change and are linked to climate anomalies (e.g., El Niño Southern Oscillation). The final shoreline position is reconstructed as the sum of oscillations at all time-scales and compares well with data from two beaches.
With recent studies forecasting potential significant change to shoreline erosion patterns, more robust prediction methods that consider changes in weather and climate patterns are increasingly critical. The approach proposed is useful because it does not depend on preexisting conditions and individual timescales can be predicted separately, accounts for nonstationarity on the drivers and shoreline response and for additional drivers that might be locally relevant.

Data Availability Statement
Data from Tairua are available at https://coastalhub.science/data. Sydney wave data was provided by the Manly Hydraulics Laboratory on behalf of the NSW Department of Planning, Industry and Environment and is available via the Australian Ocean Data Network (https://portal.aodn.org.au/). ARGUS images from which shorelines were derived are provided by the Water Research Laboratory, UNSW Australia (http://ci.wrl.unsw.edu.au/).