M. D 'anna

D Idier

B Castelle

J Rohmer

L Cagigal

F J Mendez

Maurizio D'anna
email: m.danna@brgm.fr

Effects of Stochastic Wave Forcing on Equilibrium Shoreline Modelling Across the 21 st Century Including Sea-Level Rise

Keywords: Ensemble wave forcing, shoreline modelling, equilibrium shoreline models, sea-level rise, Truc Vert

Coastal communities are currently facing the challenge of climate change and coastal retreat. While scientists are moving towards ensemble-modelling approaches to address uncertainties on shoreline evolution predictions, they rarely account for the stochastic nature of wave conditions across a variety of temporal scales (e.g., daily, weekly, seasonal, and interannual). In this contribution, we investigate the effects of the inherent variability of wave conditions on past and future multi-decadal shoreline evolution at the, cross-shore transport dominated, beach of Truc Vert (France). Using a climate-based wave emulator and variance decomposition method, we address the relative impacts of uncertain wave chronology, sea-level rise and model free parameters on modelled shoreline change, while accounting for possible correlations and interactions among the input variables. This work is done for two different wave-driven equilibrium models. The results show that the equilibrium shoreline models respond differently to the ensemble wave forcing, with strong implications on the long-term variability of modelled shoreline. We find that the modelled shoreline variance is primarily driven by the uncertain wave chronology until mid-21 st century, while the uncertainties on future sea-level rise become dominant after 2060 in all the simulated scenarios. We also found that interactions and correlations among the uncertain variables can affect the estimation of shoreline predictions uncertainties. Finally, we provide a perspective on the application of non-stationary wave-related model parameters as future research avenue for understanding uncertainties in modelled shoreline.

Introduction

Sandy beaches represent about one third of the Earth's ice-free coasts [START_REF] Luijendijk | The State of the World's Beaches[END_REF], and provide an important natural [START_REF] Schlacher | Sandy beaches at the brink[END_REF] and socio-economical resource [START_REF] Ghermandi | A global map of coastal recreation values: Results from a spatially explicit meta-analysis[END_REF] to coastal communities as well as a buffer zone during storm events. Open sandy coasts are ubiquitous worldwide, with shorelines continuously adapting to climate drivers over multiple time scales [START_REF] Montaño | A Multiscale Approach to Shoreline Prediction[END_REF][START_REF] Stive | Variability of shore and shoreline evolution[END_REF], so that temporal shifts in sea level and wave climate can strongly affect the dynamics of sandy shorelines [START_REF] Ranasinghe | Assessing climate change impacts on open sandy coasts: A review[END_REF]. On this type of coasts, shoreline change is often dominated by wave energy variations on the time scales of storms to seasons and years due to interannual variability of winter wave energy [START_REF] Dodet | Beach recovery from extreme storm activity during the 2013-14 winter along the Atlantic coast of Europe[END_REF], while sea-level rise (SLR) affects shoreline trends over longer time scales (decades to centuries) [START_REF] Ranasinghe | On the need for a new generation of coastal change models for the 21st century[END_REF]. The complex dynamics of sandy beaches together with the limited predictability of the future climate make shoreline predictions challenging, forcing coastal scientists, managers and engineers to embrace uncertainties [START_REF] Ranasinghe | On the need for a new generation of coastal change models for the 21st century[END_REF][START_REF] Toimil | Addressing the challenges of climate change risks and adaptation in coastal areas: A review[END_REF]. Indeed, the input parameters and variables of (statistical or numerical) shoreline models include uncertainties, which derive from simplifications and limited knowledge of some physical processes (model assumptions) or from the inherent variability of a model forcing (climate unpredictability), which cascade through the model and result in uncertain model predictions [START_REF] Toimil | Visualising the Uncertainty Cascade in Multi-Ensemble Probabilistic Coastal Erosion Projections[END_REF]. Thus, current practices for coastal impact assessments are moving towards probabilistic frameworks for a robust risk-informed decision-making [START_REF] Hinkel | Meeting User Needs for Sea Level Rise Information: A Decision Analysis Perspective[END_REF][START_REF] Wainwright | Moving from deterministic towards probabilistic coastal hazard and risk assessment: Development of a modelling framework and application to Narrabeen Beach[END_REF].

Probabilistic, mostly ensemble-based, approaches have become increasingly popular to study the impact several sources of uncertainties on shoreline predictions. However, most studies focused on the impact of uncertainties in future SLR on shoreline predictions while accounting for probabilistic extrapolated trends of wave-driven shoreline change [START_REF] Athanasiou | Uncertainties in projections of sandy beach erosion due to sea level rise: an analysis at the European scale[END_REF][START_REF] Cozannet | Uncertainties in Sandy Shorelines Evolution under the Bruun Rule Assumption[END_REF][START_REF] Thiéblemont | Deep uncertainties in shoreline change projections: an extra-probabilistic approach applied to sandy beaches[END_REF][START_REF] Thiéblemont | Likely and High-End Impacts of Regional Sea-Level Rise on the Shoreline Change of European Sandy Coasts Under a High Greenhouse Gas Emissions Scenario[END_REF][START_REF] Vousdoukas | Sandy coastlines under threat of erosion[END_REF], or resolving ~hourly shoreline response to deterministic realizations of future wave climate (D'Anna et al., 2021a[START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF], and did not include the inherent variability of wave climate.

The natural variability of the wave conditions introduces an unavoidable component of uncertainty to wave projections, which makes the intra-seasonal and intra-storm chronology of wave events hardly predictable [START_REF] Mankin | The Value of Initial Condition Large Ensembles to Robust Adaptation Decision-Making[END_REF]. Even when associating a comparable cumulative wave energy, different storm sequences can result in very different shoreline responses [START_REF] Baldock | Methodology for estimating return intervals for storm demand and dune recession by clustered and non-clustered morphological events[END_REF][START_REF] Coco | Beach response to a sequence of extreme storms[END_REF][START_REF] Eichentopf | Influence of storm sequencing on breaker bar and shoreline evolution in large-scale experiments[END_REF]Splinter et al., 2014a). Thus, the intrinsic uncertainty associated with the variability of wave conditions can have a significant impact on the confidence in modelled shoreline projections [START_REF] Toimil | Visualising the Uncertainty Cascade in Multi-Ensemble Probabilistic Coastal Erosion Projections[END_REF][START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF].

As the limited predictability of the future climate (and waves) introduces a certain uncertainty [START_REF] Deser | Certain Uncertainty: The Role of Internal Climate Variability in Projections of Regional Climate Change and Risk Management[END_REF] to shoreline predictions, the use of large ensembles is advocated to quantify its contribution to the uncertainty in model results [START_REF] Mankin | The Value of Initial Condition Large Ensembles to Robust Adaptation Decision-Making[END_REF]. [START_REF] Davidson | Annual prediction of shoreline erosion and subsequent recovery[END_REF] developed a method to generate ensembles of annual wave time series by sampling and shuffling monthly wave sequences from a pool of historical wave data. Although this method preserves the seasonal variability of the wave climate, it does not ensure a realistic variability on longer time scales. The recent development of stochastic climate-based wave emulators enabled the efficient generation of large ensembles of indefinitely long wave series that are characterized by different chronologies of events while obeying to realistic climate patterns [START_REF] Anderson | Time-Varying Emulator for Short and Long-Term Analysis of Coastal Flood Hazard Potential[END_REF][START_REF] Antolínez | A multiscale climate emulator for long-term morphodynamics (MUSCLE-morpho)[END_REF][START_REF] Cagigal | A multivariate, stochastic, climate-based wave emulator for shoreline change modelling[END_REF][START_REF] Pringle | A New Approach for the Stochastic Simulation of Regional Wave Climates Conditioned on Synoptic-Scale Meteorology[END_REF][START_REF] Rueda | A multivariate extreme wave and storm surge climate emulator based on weather patterns[END_REF].

Resolving long-term wave-driven shoreline change in ensemble settings requires high computational efficiency. In this context, equilibrium shoreline models are particularly convenient tools (Montaño et al., 2020) as they simulate cross-shore shoreline response to wave-driven sediment transport processes on time scales from hours to decades with low computational effort (Alvarez-Cuesta et al., 2021a;[START_REF] Davidson | Forecasting coastal evolution on time-scales of days to decades[END_REF][START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF][START_REF] Jaramillo | A shoreline evolution model considering the temporal variability of the beach profile sediment volume (sediment gain / loss)[END_REF][START_REF] Lemos | Equilibrium modeling of the beach profile on a macrotidal embayed low tide terrace beach[END_REF][START_REF] Robinet | A reduced-complexity shoreline change model combining longshore and cross-shore processes: The LX-Shore model[END_REF]Splinter et al., 2014b;[START_REF] Toimil | Managing coastal erosion under climate change at the regional scale[END_REF][START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF]. These models are based on the concept that changes in incident wave energy drive the shoreline towards a time-varying equilibrium position [START_REF] Wright | Morphodynamic variability of surf zones and beaches: A synthesis[END_REF] that is typically formulated as a function of either the current shoreline position [START_REF] Miller | A simple new shoreline change model[END_REF][START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF] or the past wave conditions [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF]Splinter et al., 2014b). Although both equilibrium approaches have shown good skill in simulating multi-year shoreline change [START_REF] Castelle | Equilibrium shoreline modelling of a high-energy meso-macrotidal multiple-barred beach[END_REF][START_REF] Montaño | Blind testing of shoreline evolution models[END_REF], their different responses to wave energy variability (D'Anna et al., 2021a;[START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF] motivate further considerations on the respective skills in reproducing shoreline change under different wave forcing scenarios, particularly on multidecadal timescales. Further, equilibrium shoreline models rely on the parametrization of some physical processes that require site-specific calibration against shoreline observations, of which the quality and availability introduce uncertainties in the model parameters [START_REF] Splinter | How much data is enough? The importance of morphological sampling interval and duration for calibration of empirical shoreline models[END_REF] that can significantly affect model predictions [START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF].

To date, few studies investigated the influence of uncertainties associated with short-and long-term variability of wave conditions on shoreline modelling. [START_REF] Cagigal | A multivariate, stochastic, climate-based wave emulator for shoreline change modelling[END_REF] applied a stochastic wave emulator to force ensemble long-term shoreline predictions to several sites dominated by cross-shore sediment transport processes to replicate shoreline return periods observed in the past. [START_REF] Kroon | Ranking uncertainty: Wave climate variability versus model uncertainty in probabilistic assessment of coastline change[END_REF] and [START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF] quantified the relative contributions of the uncertainties associated with the inherent wave variability and model parameters to longshore (1 year) and cross-shore (8 years) shoreline projections, respectively. Alvarez-Cuesta et al. (2021b) examined short-and long-term drivers of the mean shoreline change at two longshore transport dominated sites using a wave ensemble and three percentiles of future SLR for two RCP scenarios. [START_REF] Toimil | Visualising the Uncertainty Cascade in Multi-Ensemble Probabilistic Coastal Erosion Projections[END_REF] produced multi-ensemble 2081-2100 projections to quantify the uncertainties associated with the individual steps of the modelling process, including uncertainties on SLR, waves, storm surges, and future climate scenario for two shoreline models. The works mentioned above resolved cross-shore wave-driven shoreline change using disequilibrium approach of the same kind, that is, where the beach equilibrium state is defined by the current shoreline position [START_REF] Miller | A simple new shoreline change model[END_REF][START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF].

Many efforts have been dedicated to quantify the impact of intrinsic uncertainty on shoreline predictions.

However, the role of the inherent wave conditions variability on ensemble shoreline predictions in relation to uncertain SLR and model free parameters, and how substantially different modelling approaches respond to this type of uncertainties, has not been investigated. In addition, the wave-driven model free parameters uncertainties are strongly connected to the wave forcing variability (D'Anna et al., 2021a[START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF][START_REF] Davidson | Annual prediction of shoreline erosion and subsequent recovery[END_REF][START_REF] Ibaceta | Enhanced Coastal Shoreline Modeling Using an Ensemble Kalman Filter to Include Nonstationarity in Future Wave Climates[END_REF][START_REF] Montaño | A Multiscale Approach to Shoreline Prediction[END_REF][START_REF] Splinter | Rapid adjustment of shoreline behavior to changing seasonality of storms: observations and modelling at an open-coast beach[END_REF]Splinter et al., , 2014b;;[START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF], and the effects of possible interactions and correlations between the stochastic wave forcing and other uncertain model inputs have not yet been investigated.

In this contribution, we investigate the response of two different wave-driven equilibrium shoreline models to wave-forcing ensembles, and investigate how stochastic wave forcing affects the respective long-term shoreline predictions in relation to other uncertain factors (e.g. SLR) while accounting for possible interactions among them. A Global Sensitivity Analysis (GSA, [START_REF] Saltelli | Global sensitivity analysis: The premier[END_REF] is performed to quantify the contributions of uncertain variables to the total uncertainties on shoreline trajectories modelled with two equilibrium modelling approaches. The GSA estimates the model sensitivity 'globally' by exploring all the range of plausible values for the uncertain inputs while accounting for all possible interactions and statistical dependence between them. The analysis is carried out on the cross-shore transport dominated Truc Vert beach (France) accounting for uncertainties in free model parameters, future wave climate variability and SLR projections for two future representative greenhouse-gas concentration pathway (RCP) scenarios. We investigate an additional scenario over the past 23.5 years, where only wave chronology and model free parameters are considered as uncertain input variables.

The remainder of the paper provides a description of Truc Vert beach, the data, the shoreline models and the method used herein (Section 2); the assessment of the probability distributions associated to the uncertain model inputs, required for the GSA, and the ensemble modelling setup (Section 3). The results are presented and discussed in Section 4 and 5, respectively, followed by the conclusions (Section 6).

Study Site, Material and Method

Truc Vert beach

Truc Vert is a sandy beach located within the ~130 km long straight Gironde coast, southwest France.

It is backed by a high (~20 m) and large (~250 m) aeolian dune [START_REF] Robin | 150 years of foredune initiation and evolution driven by human and natural processes[END_REF] separating the beach from a large area of state-owned forest (Figure 1a-d). The beach is exposed to a highly energetic wave climate generated in the North Atlantic Ocean. Incident wave energy mostly comes from the westnorthwest direction with a strong seasonal modulation of significant wave height (Hs) and peak wave period (Tp) (Castelle et al., 2018a;[START_REF] Charles | Climate change impact on waves in the Bay of Biscay, France[END_REF][START_REF] Robinet | Statistical modeling of interannual shoreline change driven by North Atlantic climate variability spanning 2000-2014 in the Bay of Biscay[END_REF] reaching monthly mean values of 1.1 m and 9 s in August, and 2.4 m and 12.8 s in January, respectively. The local winter wave energy also shows interannual variations driven by natural large-scale patterns of atmospheric variability such as the West Europe Pressure Anomaly [START_REF] Castelle | A new climate index controlling winter wave activity along the Atlantic coast of Europe: The West Europe Pressure Anomaly[END_REF], which affects the spatio-temporal distribution of storm events and clustering winters [START_REF] Castelle | Impact of the winter 2013-2014 series of severe Western Europe storms on a double-barred sandy coast: Beach and dune erosion and megacusp embayments[END_REF][START_REF] Masselink | Extreme wave activity during 2013/2014 winter and morphological impacts along the Atlantic coast of Europe[END_REF]. Crossshore processes driven by the temporal variability of incident wave energy have been recognized to be the main driver of shoreline changes at Truc Vert [START_REF] Castelle | Equilibrium shoreline modelling of a high-energy meso-macrotidal multiple-barred beach[END_REF][START_REF] Robinet | A reduced-complexity shoreline change model combining longshore and cross-shore processes: The LX-Shore model[END_REF][START_REF] Robinet | Statistical modeling of interannual shoreline change driven by North Atlantic climate variability spanning 2000-2014 in the Bay of Biscay[END_REF], in line with the negligible longshore drift gradients in this area [START_REF] Idier | Longshore sediment flux hindcast: spatiotemporal variability along the SW Atlantic coast of France[END_REF]. Truc Vert is a mesomacro tidal beach with a ~3.7 m annual mean spring tidal range and a maximum astronomical tide range reaching up to 5 m (Castelle et al., 2018a), with negligible tidal currents compared to wave-driven currents in the nearshore zone.

Starting from the early 2000s, Truc Vert beach has been monitored monthly to bi-monthly with topographic DGPS surveys from 2003 (with 1-year gap in 2008) (Figure 1e), with additional field campaigns including a 5-weeks campaign of daily bathymetric surveys in 2008 [START_REF] Parisot | Treatment of topographic and bathymetric data acquired at the Truc-Vert Beach during the ECORS Field Experiment[END_REF][START_REF] Sénéchal | The ECORS-Truc Vert'08 nearshore field experiment: presentation of a three-dimensional morphologic system in a macro-tidal environment during consecutive extreme storm conditions[END_REF], and seasonal high-resolution digital elevation models from photogrammetry of UAV images covering 4 km of beach-dune since 2017 [START_REF] Laporte-Fauret | Low-Cost UAV for High-Resolution and Large-Scale Coastal Dune Change Monitoring Using Photogrammetry[END_REF]. Here, following [START_REF] Castelle | Equilibrium shoreline modelling of a high-energy meso-macrotidal multiple-barred beach[END_REF] and [START_REF] Robinet | Statistical modeling of interannual shoreline change driven by North Atlantic climate variability spanning 2000-2014 in the Bay of Biscay[END_REF], the shoreline position is defined as the intersection of the beach profile with the 1.5 m level contour above mean sea level, which approximately corresponds to the mean high water level at Truc Vert beach and best correlates to variations in total beach-dune volume. For more details on the topo-bathymetric datasets, see [START_REF] Castelle | 16 years of topographic surveys of ripchannelled high-energy meso-macrotidal sandy beach[END_REF].

Over the past few decades, the Truc Vert coastline as well as the dune foot have been observed to be approximately stable (Castelle et al., 2018b), with some exceptions such as the extreme winter of 2013-2014, characterized by rapid successions of severe storms that caused large erosion and dune scarping [START_REF] Castelle | Impact of the winter 2013-2014 series of severe Western Europe storms on a double-barred sandy coast: Beach and dune erosion and megacusp embayments[END_REF][START_REF] Masselink | Extreme wave activity during 2013/2014 winter and morphological impacts along the Atlantic coast of Europe[END_REF]. The absence of coastal or inland structures nearby, together with the wealth of available data, make this site a worthy benchmark for modelling and interpreting shoreline response to natural processes.

Sea level and wave data

Historical wave data was required for calibration of the shoreline models and for generating the wave ensemble used in this work to force the probabilistic past and future shoreline projections. Information on past and future relative mean sea level (MSL) at Truc Vert beach was needed to account for SLRdriven shoreline erosion during model calibration and simulations. Herein we used the same datasets of past wave conditions, past MSL, and future SLR projections adopted in D' Anna et al. (2021a). In this Section, we briefly describe these datasets, which are all extended here to June 2021.

Historical incident wave conditions at Truc Vert beach were obtained from the NORGAS-UG regional wave model [START_REF] Michaud | Improvements of the new French coastal wave forecasting system and application to a wave-current interaction study[END_REF]. The NORGAS-UG model provides hourly hindcast wave conditions over the French Atlantic coastal area, and has been extensively validated against data from a number of wave buoys [START_REF] Michaud | Improvements of the new French coastal wave forecasting system and application to a wave-current interaction study[END_REF]. We extracted the time series of significant wave height (Hs), peak period (Tp) and mean direction (Dm) from January 1994 to June 2021, at the model grid point co-located with the in situ CANDHIS directional wave buoy (44°39′9″N; -1°26′48″W) moored in 54-m depth (Figure 1a). Figure 2a shows the characteristic strong seasonality of the incident wave energy (synthetized here by H²sTp) at Truc Vert, with also prominent winter wave energy interannual variability. We adopted a constant rate of historical relative SLR of 3.3 ± 0.7 mm/year (median ± σ), obtained by D' Anna et al. (2021a) from the combination of geocentric MSL change and vertical land motion at Truc Vert beach. Future relative MSL at Truc Vert beach was estimated using global SLR projections until 2100 and the respective likely ranges (17 th -83 rd percentiles) from the Special Report on Oceanic and Cryosphere in a Climate Change (SROCC; [START_REF] Oppenheimer | Sea Level Rise and Implications for Low-Lying Islands, Coasts and Communities[END_REF], for the RCP4.5 and RCP8.5 scenarios. The global MSL estimates are projected to Truc Vert beach accounting for the regional fingerprint of each process contributing to SLR, including Glacial Isostatic Adjustment [START_REF] Slangen | Projecting twenty-first century regional sea-level changes[END_REF], in line with [START_REF] Thiéblemont | Likely and High-End Impacts of Regional Sea-Level Rise on the Shoreline Change of European Sandy Coasts Under a High Greenhouse Gas Emissions Scenario[END_REF]. Finally, the contribution of vertical land motion was added assuming that this process will maintain the past observed rate of 1.2 ± 0.6 mm/yr over the 21 st century (as in D'Anna et al., 2021a).

Shoreline change models

As Truc Vert shoreline change is mostly driven by cross-shore sediment transport processes (Castelle et al., 2014;[START_REF] Robinet | A reduced-complexity shoreline change model combining longshore and cross-shore processes: The LX-Shore model[END_REF][START_REF] Robinet | Statistical modeling of interannual shoreline change driven by North Atlantic climate variability spanning 2000-2014 in the Bay of Biscay[END_REF], we resolved cross-shore shoreline response to the incident wave climate using the ShoreFor [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF]Splinter et al., 2014b) and the [START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF] equilibrium shoreline models. Given that the Truc Vert offshore bathymetric iso-contours are mostly shore-parallel, breaking wave conditions, which are required to force the above models, were estimated with the [START_REF] Larson | Direct Formula to Compute Wave Height and Angle at Incipient Breaking[END_REF] empirical formula. Long-term shoreline retreat induced by SLR was modelled using the [START_REF] Bruun | Sea-level rise as a cause of shore erosion[END_REF] Rule. The wave-driven models and the Bruun Rule were run separately and then linearly combined in order to avoid spurious feedback mechanisms between the models (D'Anna et al., 2021b), in line with previous works [START_REF] Vitousek | A model integrating longshore and cross-shore processes for predicting long-term shoreline response to climate change[END_REF].

Wave-driven shoreline models, disequilibrium conditions and free parameters

The ShoreFor and [START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF] models are both rooted in the equilibrium beach concept, which assumes the existence of a time-dependent equilibrium state between the beach and the incident wave conditions. Based on this concept, the rate and direction of shoreline change depends on the combination of the incident wave force and the disequilibrium state of the beach [START_REF] Wright | Morphodynamic variability of surf zones and beaches: A synthesis[END_REF]. The main difference between ShoreFor and the [START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF] models lies in the approach adopted to quantify the (dis)equilibrium state and how this responds to the time-variability of wave conditions.

ShoreFor

The ShoreFor model (herein SF) disequilibrium approach is based on the present and past dimensionless fall velocities (Ω), and the concept of beach memory. The rate of shoreline change is estimated as: 𝑑𝑌 𝑑𝑡 = 𝑘 𝑠 +/-𝑃 0.5 ∆𝛺 𝜎 ∆𝛺 (1) where, ks +/-(m s -1 W -0.5 ) is a response rate parameter, P(W) is the wave power at breaking, ΔΩ/σΔΩ is the disequilibrium term, and b(m/s) is a linear trend term. The disequilibrium state at a given time is defined as the normalized difference between the equilibrium and current dimensionless fall velocities (ΔΩ = Ωeq -Ω), here calculated from offshore wave conditions [START_REF] Robinet | A reduced-complexity shoreline change model combining longshore and cross-shore processes: The LX-Shore model[END_REF]:

𝛺 = 𝐻 𝑠,𝑜 𝑇 𝑝 𝑤 (2) 𝛺 𝑒𝑞 = ∑ 𝛺 𝑖 10 -1 𝜑 ⁄ 2𝜑 𝑖=1 [∑ 10 -1 𝜑 ⁄ 2𝜑 𝑖=1 ] -1 (3)
where w(m/s) is the sediment fall velocity, and φ(days) is a beach memory parameter. The equilibrium term (Ωeq, Equation 3) is defined by the weighted average of Ω over the past 2φ days, attributing to φ the role of 'memory' of the beach. As w is generally assumed constant over time and φ is a stationary model parameter, the dynamic equilibrium state fully relies on the variability of wave conditions, providing SF with the ability to respond to long-term variations of wave regime.

The response rate parameter ks +/-represents an efficiency rate, which determines the rate of shoreline change in response to the wave forcing. This parameter associates different values for accretion (ΔΩ>0, ks + ) and erosion (ΔΩ<0, ks -) events to account for the different time scales related to the respective driving processes. In SF, the ks values for accretion and erosion are assumed to be linearly related by a so-called 'erosion rate' parameter r (ks -= r ks + ). The r parameter is such that an increasing (decreasing) trend the wave forcing over the simulated period produces an erosive (accreting) trend of shoreline change, and is estimated as follows:

𝑟 = | ∑ 〈𝐹 + 〉 𝑁 𝑖=1 ∑ 〈𝐹 -〉 𝑁 𝑖=1 | (4) 𝐹 = 𝑃 0.5 ∆𝛺 𝜎 ∆𝛺 ( 5
)
where F + and F -refer to accretion (ΔΩ >0) and erosion (ΔΩ <0) events, respectively, 〈.〉 indicates an operation that removes the linear trend, and N is the number of time-steps over the simulated period.

The parameter b represents a residual, steady, shoreline trend due to any unresolved long-term processes over the simulated period (e.g. SLR, longshore sediment transport gradients, etc.). φ, ks + and b are sitespecific free parameters that require calibration against observed wave and shoreline data. Given the negligible contribution of long-term processes such as longshore drift gradients in the Truc Vert area [START_REF] Idier | Longshore sediment flux hindcast: spatiotemporal variability along the SW Atlantic coast of France[END_REF], and that uncertainties on the linear trend explode over the long time scales [START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF] simulated here (decades to centuries), we set b=0.

Yates Model

In the [START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF] , ΔE < 0). However, contrary to SF, in Y09 ky + , and ky -are independent free parameters. The equilibrium wave energy Eeq is related to the shoreline position through the empirical equilibrium energy function:

𝐸 𝑒𝑞 (𝑌) = 𝑎 1 𝑌 + 𝑎 2 (7)
where a1(m 2 /m) and a2(m 2 ) are empirical coefficients. While a1 and a2 associate the dimensions of 'energy per meter' and 'energy', respectively, these are to be interpreted as empirical parameters. In Y09, the values for a1, a2, ky + , and ky -are specific to the site of application, and need calibration against available wave and shoreline observations. Y09 does not include an explicit 'beach memory' term, but implicitly accounts for recent past events through Equation 7. In fact, this relationship allows damping the efficiency of wave-driven shoreline change as the shoreline approaches its equilibrium position [START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF][START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF]. However, as the coefficients of Equation 7 (a1, a2) are estimated based on the wave conditions observed during the calibration period and remain stationary over time, Equation 7does not respond to shifts in wave climate. Recently, [START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF] showed that a rearrangement of Y09 free parameters results in combined terms that can be interpreted as local equilibrium time and spatial scale factors. However, in this application we retained the original form of Y09, interpreting a1 and a2 as empirical parameters.

Sea-level induced shoreline erosion

The chronic shoreline recession driven by SLR was estimated here with the [START_REF] Bruun | Sea-level rise as a cause of shore erosion[END_REF] model. This model is based on cross-shore sediment balance conservation under the assumption that the equilibrium beach profile shape remains unchanged over time. According to the Bruun model, over time scales larger than years, the beach profile adapts to the slow-rising sea level with an upward and landward shift resulting in a mean shoreline retreat (dYSLR) given by:

𝑑𝑌 𝑆𝐿𝑅 𝑑𝑡 = 𝑆𝐿𝑅 𝑟𝑎𝑡𝑒 𝑡𝑎𝑛𝛽 (8)
where SLRrate is the rate of SLR(m/time) and tanβ is the mean slope of the active beach profile. Here, the latter profile was defined between the dune crest and the 'depth of closure', beyond which no significant sediment exchange is assumed to occur [START_REF] Bruun | The Bruun Rule of Erosion by Sea-Level Rise: A Discussion on Large-Scale Twoand Three-Dimensional Usages[END_REF][START_REF] Wolinsky | A unifying framework for shoreline migration: 2. Application to wave-dominated coasts[END_REF]. We estimated tanβ = 0.023 using the beach profile reported in Figure 1c and the depth of closure calculated according to [START_REF] Hallermeier | Uses for a calculated limit depth to beach erosion[END_REF]. Although the generalized applicability of the Bruun model has been questioned due to its restricting assumptions [START_REF] Cooper | Sandy beaches can survive sea-level rise[END_REF][START_REF] Cooper | Sea-level rise and shoreline retreat: time to abandon the Bruun Rule[END_REF][START_REF] Ranasinghe | Estimating coastal recession due to sea level rise: beyond the Bruun rule[END_REF], Truc Vert beach settings are in line with most of Bruun's underlying assumptions (D'Anna et al., 2021b).

Method

Multivariate stochastic climate-based wave emulator

Equilibrium shoreline models require continuous time series of incident wave conditions with a time resolution of the order of hours. In this study, we produced an ensemble of synthetic wave forcing conditions using the climate-based stochastic wave emulator developed by [START_REF] Cagigal | A multivariate, stochastic, climate-based wave emulator for shoreline change modelling[END_REF]. For a given location, the emulator exploits broad spatial scale sea level pressure fields and synoptic patterns, together with local historical wave data, to generate ensembles of unique continuous realistic wave time series. The wave series are generated through: (i) a statistical downscaling of climate forcing relying on annual and daily synoptic patterns; (ii) an intra-seasonal chronology model based on autoregressive logistic regression [START_REF] Anderson | Time-Varying Emulator for Short and Long-Term Analysis of Coastal Flood Hazard Potential[END_REF]; and (iii) a 'stretching and shuffling' of intra-storm historical wave properties using Gaussian copulas. This method ensures that the modelled wave time series present different chronologies of events while preserving intra-storm to interannual statistical properties consistent with the historical data. The latter means that the generated time series are also characterized by the same long-term trend observed over the historical wave data. The wave emulator has already been applied to produce ensembles of wave time series for probabilistic shoreline modelling at several sites based on the respective records of historical wave conditions, and showed good skill in reproducing observed shoreline erosion statistics [START_REF] Cagigal | A multivariate, stochastic, climate-based wave emulator for shoreline change modelling[END_REF][START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF]. For a more detailed outline of the underlying methodology see [START_REF] Cagigal | A multivariate, stochastic, climate-based wave emulator for shoreline change modelling[END_REF].

We used the wave emulator to model 200 hourly time series of wave conditions (Hs, Tp, Dm) from January 1994 to December 2099, based on 27 years of wave hindcast data offshore Truc Vert beach (Section 2.2) and sea level pressure fields from the Climate Forecast System Reanalysis [START_REF] Saha | The NCEP Climate Forecast System Reanalysis[END_REF]. Figure 2b shows the mean (orange line), standard deviation (yellow area), and the min-max envelope (red area) of the emulated Hs series, which attained extreme values of same order of magnitude as observed during the hindcast period (Hs ~8-10 m). Figure 2c illustrates the standard deviation (yellow area) and envelope (red area) of the 3-month averaged Hs 2 Tp ensemble, along with an example of individual time series (black line). The statistical properties of the generated ensemble showed good agreement with the historical wave data. Figure 2d,e show a quantile-quantile comparison of Hs between the reference hindcast data and each member of the wave ensemble for the full time series and their seasonal repartition, respectively. Here, seasons were defined as follows: Winter (December-January-February), Spring (March-April-May), Summer (June-July-August) and Autumn (September-October-November). The comparisons of Tp and Dm are shown in the Supplementary Material (Figure S1). While the Hs quantiles estimated on the entire time series show a near-perfect consistency (Figure 2c), the seasonal plots show some discrepancies (Figure 2d). However, these differences are not expected to affect the analyses carried out in this work (as discussed in Section 5.4). Consistently with the historical wave series, the winter mean wave heights of the emulated time series show no significant trend between present day (2021) and the end of the 21 st century (Figure S2).

Global Sensitivity Analysis

The uncertainties associated to model assumptions and climate unpredictability propagate through the model resulting in uncertain model predictions. Regardless of their nature, the uncertainties affecting the inputs can have a different influence on the model outcome. In variance-based approaches, the uncertainties associated with a variable are represented by its variance. We quantified the different contributions to the uncertainties on modelled shoreline change at Truc Vert beach by performing a variance-based GSA [START_REF] Saltelli | Global sensitivity analysis: The premier[END_REF]. The variance-based GSA explores the full range of realistic combinations of uncertain model inputs with the aim to divide the uncertainties (i.e. variance) of probabilistic model results into several portions, each one attributed to an uncertain model input (or a set of inputs). The resulting contribution that an uncertain input Xi provides to the uncertainties in model predictions is given by the main effect (first-order sensitivity), and by interactions terms (higher-order sensitivity). The main effect measures the sensitivity of model results to variations of the Xi value alone, while the higher-order components measure the additional effects that simultaneous variations of Xi with other uncertain inputs have on the results' variance. The latter accounts for the fact that simultaneous variation of multiple Xi's may drive a larger contribution to the results variance than the cumulated main effects (individual variations) of Xi's. For each uncertain input Xi, the normalized main effect of Xi uncertainties on the model results can be synthetized by a sensitivity measure (S'i), also known as firstorder Sobol' indices [START_REF] Sobol′ | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], expressed as:

𝑆′ 𝑖 = Var(E(𝑌|𝑋 𝑖 )) Var(𝑌) (9)
where Y is the model prediction, Xi is the i-th uncertain input variable, E(Y|Xi) is the expected value of Y conditional to Xi, and Var(E(Y|Xi)) is the variance of the E(Y|Xi)s obtained for all the possible values of Xi. S'i quantifies the percentage of the model results variance that can be attributed to the uncertainties on the variable Xi alone. When Xi interacts with other uncertain inputs (e.g. Xj) within the model, the uncertainties associated with Xj may amplify the impact of Xi on the results uncertainties (and those of Xi may amplify the impact of Xj) compared to the first-order effects. In this case, the contribution of Xi and Xj to the variance of model results is larger than S'i+S'j. The latter contribution as well as higherorder effects can be estimated by calculating higher-order Sobol' indices with formulations analogous to Equation 9 (see Text S2 of Supplementary Material).

The GSA described above builds on the assumption that all the uncertain input variables are statistically independent from one to another. If two or more uncertain inputs are statistically dependent, say Xi and Xj, such correlations can significantly influence the uncertainties of model results, and the S'i includes some information that is also included in S'j, regardless of the possible interactions between the two variables [START_REF] Do | Correlation Effects? A Major but Often Neglected Component in Sensitivity and Uncertainty Analysis[END_REF][START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with sobol' indices, numerical estimation and applications[END_REF]. This is the case, for instance, of the calibrated shoreline model parameters, which can associate up to more than 70% correlation coefficient (see Section 3.1).

The Shapley Effects (Sh's) have been proposed based on the Shapley values [START_REF] Shapley | A Value for n-Person Games[END_REF] to overcome the difficulties in the interpretation and use of the Sobol' indices in presence of statistical dependence and interaction among inputs [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: Comparisons with sobol' indices, numerical estimation and applications[END_REF][START_REF] Owen | Sobol' Indices and Shapley Value[END_REF][START_REF] Song | Shapley Effects for Global Sensitivity Analysis: Theory and Computation[END_REF]. The Shapley values were originally introduced in game theory to evaluate the "fair share" of team players after a game based on their individual contributions and their interactions with other players. In the context of the GSA the 'players' are the uncertain model inputs, and the model results variance is the 'value' to share. The Shi quantifies the percentage of output variance (0≤ Shi ≤1) associated to the input variable Xi, including the main effect, an equitable share of all its interaction terms with other variables, and the contributions of possible statistical dependence (correlation) among the uncertain inputs, and is expressed as follows: Sh's include the first-and higher-order terms and the shared effects of inputs correlations to the model results, they do not provide these contributions separately.

𝑆ℎ 𝑖 = 1 𝑘 ∑ ( 𝑘 -1 |𝐴| ) -1
Here, we analysed the time-varying effects of uncertainties associated to model free parameters, SLR, and wave chronology on the modelled Truc Vert shoreline change using the Shapley Effects, while firstorder Sobol' indices were computed to support the discussion in Section 5.2. In addition, we assessed the impact of the shoreline modelling approach choice by repeating the application for each of the two wave-driven models (SF and Y09). The framework described above requires the definition of a probability distribution for each uncertain input variable, a large number of realizations based on a Monte-Carlo procedure (accounting for statistical dependence among the variables), and the computation of Shapley effects at each analysed time step (Figure 3).

Input Uncertainties

Model free parameters

Equilibrium shoreline models are affected by uncertainties related to the underlying physical assumptions as well as their data-driven nature. We accounted for the former by running and comparing the GSA for different scenarios using the disequilibrium approaches of the SF and Y09 models (Section 2.3.1). These models require the calibration of free parameters that associate uncertainties depending on the extent, density and quality of the available shoreline datasets [START_REF] Splinter | How much data is enough? The importance of morphological sampling interval and duration for calibration of empirical shoreline models[END_REF]. The uncertainties on the calibration of SF (ks + , φ) and Y09 (ky +/-, a1, a2) model free parameters were represented by the respective joined probability distributions, defined using the approach adopted by D' Anna et al. (2021a).

In this approach, the shoreline models were calibrated using the Simulated Annealing algorithm [START_REF] Bertsimas | Simulated Annealing[END_REF], which generates a large number of model parameters combinations before converging to a set of best-fit parameters. The parameters combinations that produce a model performance higher than a predefined threshold were selected and used to fit an empirical multivariate probability distribution (multivariate kernel function). The model performance was measured in terms of Nash-Sutcliffe score (NS) [START_REF] Nash | River flow forecasting through conceptual models part I -A discussion of principles[END_REF], which quantifies the model skill compared to the skill associated with the mean of the observed shoreline position (𝑌 𝑜 ̅ ).

𝑁𝑆 = 1 - ∑ (𝑌 𝑚 𝑛 -𝑌 𝑜 𝑛 ) 2 𝑁 𝑛=1 ∑ (𝑌 𝑜 ̅ -𝑌 𝑜 𝑛 ) 2 𝑁 𝑛=1 ( 12
)
where N is the number of observations, Ym n and Yo n are the n-th modelled and observed shoreline positions, respectively. An NS value of 1 indicates that the model perfectly reproduces the observed shoreline positions, while NS=0 corresponds to a model as skilful as the straight mean (𝑌 0 ̅ ). The minimum threshold used to fit the free parameters probability distribution is set to NS ≥0.25, which is an appropriate value for SF and Y09 applications at Truc Vert beach (D'Anna et al., 2021a). Here, we calibrated the models over the period March 2012 -September 2019, which show no long-term shoreline trend supporting the assumption that b=0 in the SF model. Further available shoreline observations up to June 2021 (Figure 1e) were exploited for the validation of the ensemble simulations. Shoreline data prior to 2012 is discarded from the calibration procedure due to the lower confidence deriving from the limited alongshore coverage of their associated topographic survey [START_REF] Castelle | 16 years of topographic surveys of ripchannelled high-energy meso-macrotidal sandy beach[END_REF][START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF]. Over the calibration period, the Y09 and SF models produce a Root-Mean-Square-Error (RMSE) of 5.1 m and 7.2 m, respectively, and a coefficient of determination R 2 of 0.80 and 0.65, respectively.

The resulting probability distributions of Y09 and SF model free parameters are shown in Figure 4a,b together with the respective Pearson correlation coefficients (denoted ) indicating several levels of correlation among the model parameters. The ranges of possible parameters values are reported in Table 1.

Sea-Level Rise

Information on past relative MSL associates uncertainties related with measurement errors of the geocentric MSL and the vertical land motion. This kind of uncertainties are very small compared to the uncertainties on the free parameters of shoreline models, and have been shown to have a negligible impact on the uncertainties of shoreline change at Truc Vert over the past 20 years [START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF].

Therefore, here the uncertainties on past SLR were not included in the GSA.

Uncertainties on future SLR estimates essentially derive from unknowns on the future climate scenario and modelling of the different contributing processes. These sources of uncertainty were accounted through the likely range of the SROCC projections.

Waves chronology

The stochastic nature of wave conditions provides a portion of uncertainties to the modelled shoreline change, which responds differently to different chronologies of wave events. In order to account for the uncertainties on wave chronology, we define an indicator variable (i.e. categorical variable) that takes up discrete values between 1 and 200, each value being associated to a given time series described in Section 2.4.1 (index 1 for time series N°1, index 2 for time series N°2, etc.). Then, we randomly and uniformly sample the values of the indicator variable to select the time series used to force the shoreline modelling realizations within the GSA. In the following, the influence of the wave chronology is analysed by assessing the Shapley Effect of the defined indicator variable. This differs from the analysis of uncertain model free parameters and future SLR percentiles in that the indicator variable is not a quantitative value but a categorical number. This approach, has been applied by [START_REF] Rohmer | Combining meta-modeling and categorical indicators for global sensitivity analysis of long-running flow simulators with spatially dependent inputs[END_REF] for spatially varying uncertain inputs, allows accounting for the bulk effect of stochastic time-varying wave forcing on the uncertainties of modelled shoreline positions.

Models setup

Ensembles of simulated shoreline time series were generated for one 'past scenario' and four 'future scenarios' at Truc Vert beach. All scenarios covered the calibration period (March 2012 -September 2019) and a validation period (September 2019 -June 2021). In the 'past scenario', shoreline evolution was simulated from January 1998 to June 2021 using SF and the Bruun model. As the Y09 computation of shoreline change at a given time is dependent on the previous shoreline position, the model does not allow propagating the uncertainties on modelled shoreline towards the past. Therefore, the Y09 model application to the 'past scenario' was not included in the main analysis, as discussed in Section 5.1.

Although wave data is available on the 'past scenario' simulated period, here we considered random wave conditions characterized by the same statistical properties as the hindcast wave data, i.e. wave data generated with the stochastic climate-based wave emulator (Section 2.4.1). In the 'future scenarios', the shoreline position was projected over March 2012 -December 2099 using Y09 and SF with the Bruun model for the RCP4.5 and 8.5 SLR scenarios. The models were forced using MSL time series corresponding to percentiles randomly sampled between the 5 th and 95 th , and the randomly selected wave time series characterized by the same properties as the past wave climate (Section 3.3).

For each ensemble, several thousands of shoreline trajectories were simulated using different combinations of model free parameters, sea level percentiles and wave time series sampled from the respective probability distributions. One exception is the 'past scenario', where local SLR is considered deterministic with a constant rate of 3.31 mm/year, while only the model free parameters and forcing wave series were sampled for each model run. Over the calibration and validation periods, the models were forced using the observed (deterministic) wave and sea-level time series. The starting time of both past and future shoreline simulations was set to the starting date of the models calibration (23 rd May 2012), so that the modelled shoreline trajectories fit the shoreline data as observed during the calibration and the uncertainties on simulated past shoreline change propagate towards the past. All simulations are run with a 3-hour time step and outputs were stored with a 2-weeks frequency. The ensembles of model results were analysed through the envelope (min-max) and the likely range, defined here as the standard deviation interval (µ±σ), of the simulated shoreline positions. The former is an indicator of the shoreline response to extreme wave events, while the latter is proportional to the results variance (σ²) used to perform the GSA. Herein, positive (negative) variations of the shoreline position dY/dt>0 (dY/dt<0) indicate shoreline accretion (erosion).

The Shapley Effects (Section 2.4.2) were computed for each stored model output excluding the calibration and validation periods (i.e. excluding March 2012 -June 2021) by post-processing the inputs and the results of the ensemble model realizations with a nearest-neighbour search technique [START_REF] Broto | Variance Reduction for Estimation of Shapley Effects and Adaptation to Unknown Input Distribution[END_REF], using the 15 nearest neighbours. Before the Sh's computation, the uncertain input values were rescaled with a whitening procedure in order to achieve a homogeneous scale, and the wave time series were treated as a categorical input, i.e. they are identified as the i-th, … , j-th, time series. For each GSA application, we estimated 90% confidence bounds applying a 100-iteration bootstrapping at three different times of the simulated period for the 'future scenarios' (1 st Jan. 2030, 2060 and 2090) and two for the 'past scenario' (1 st Jan. 2000 and 2010). Confidence intervals were estimated by sampling without replacement with a proportion of the total sample size of 0.80, drawn uniformly (following the implementation in R sensitivity package).

Results

Past and future shoreline ensembles

The modelling setup described in Section 3.4 (Table 2) produced five ensembles of simulated shoreline position time series: one covering the past 23.5 years (Figure 5a), and four spanning from present day to 2100 (Figure 5b-e). 2a), with shoreline data falling within the range of the modelled ensembles (Figure 5). Compared to the 'future scenarios', the SF+B application to the past scenario shows slightly larger ranges of modelled shoreline over the calibration and validation periods (Figure 5a). This is due to the influence of the ensemble forcing wave conditions prior to 2012 on the SF shoreline change computed over the following period.

In the 'past scenario' (Figure 5a), the likely range (envelope) of shoreline position varies from [-29 m; to 2012, the Sh values of the SF parameters ks + and φ first decrease from 15% to 5% and 13% to 5% (~2006), respectively, then increase to 40% and to 39%, respectively by 2012. The wave chronology's Sh gradually increases from 78% to 85% (~2006) before decaying to 21% (Figure 6b-d). The 90% confidence bounds were estimated to be within ±1% for all the variables. All Sh's show a mild seasonal signal, where the Sh's of free parameters increase during accretion (low wave energy) periods, with a similar but out-of-phase signal for wave chronology's Sh. Figure 7 and Figure 8 show the GSA results for the 'future scenarios' where wave-driven shoreline change is simulated using the SF and Y09 models, respectively. For both modelling approaches, the Sh's show a similar behaviour between the RCP4.5 and RCP8.5 scenarios, with notable differences only in the second half of the 21 st century, where the two SLR uncertainties diverge significantly.

In the SF scenarios, the model free parameters ks + and φ control the initial model variance of the projection period, with Sh's attaining 32% and 48%, respectively, before exponentially decaying to 6% and 5% by 2030 (Figure 7b,c). Starting from 2030, both Sh's of ks + and φ gradually increase up to 10% by 2095, and further up to ~18% by 2100. The SLR Sh rapidly drops from 20% to 5% between 2021 and 2025, then grows gently back to ~20% by 2080 and dramatically up to 80% by 2100 in the RCP8.5 scenario (Figure 7d). In the RCP4.5 scenario, the model sensitivity to SLR remains below 10% until 2097 before increasing up to 23% by 2100 (Figure 7d). The impact of the uncertain wave chronology raises up to ~83% from 2021 to 2025, and then decreases progressively mirroring the behaviour of SLR Sh (Figure 7). In the RCP8.5 (RCP4.5) scenario, the wave chronology Sh decays down to ~70% (~80%) by 2080, and accelerates down to 0% (~30%) at the end of the simulated period. When using the Y09 approach, the GSA produces the series Sh shown in Figure 8. The Sh's of the empirical Y09 parameters a1 and a2 decrease regularly from the initial values of 23% and 14%, respectively, to 10% by 2100 (Figure 8b,c). Over the projection period, the Sh of ky + (ky -) decreases from 23% to 10% (10% to 8%), and associates a well-defined (weak) seasonal signal with higher values during the summer (winter) (Figure 8d,e). SLR's Sh increases regularly from ~5% in 2021 to 20% in 2050, and up to 58% and 32% in the RCP8.5 and RCP4.5 scenarios, respectively (Figure 8f). The wave chronology's Sh shows a decreasing long-term trend mirroring the behaviour of SLR's Sh and a marked seasonality. In both RCP scenarios, the waves' Sh ranges between 45% and 18% over the first decade (2021)(2022)(2023)(2024)(2025)(2026)(2027)(2028)(2029)(2030), and between 17% and 1% (23% and 1%) over 2090-2100 for the RCP8.5 (RCP4.5) scenario (Figure 8e).

Discussion

Effects of stochastic wave chronology

Past scenario

Despite the availability of wave hindcast data, we forced the 'past scenario' simulations using ensemble (stochastic) wave conditions. This choice is motivated by two reasons: (i) including stochastic wave chronologies in the 'past scenario' allows to analyse the relative impact of the intrinsic uncertainties of waves and the epistemic uncertainties of model free parameters on simulated shoreline; and (ii) since the wave forcing ensemble is based on the aforementioned hindcast data, existing shoreline hindcasts and observations can provide a qualitative validation for the emulated wave time series.

In the 'past scenario', where shoreline change is simulated using SF and the Bruun model and uncertainties on SLR are omitted, the modelled shoreline shows an overall accreting trend from 1998 to 2012. This is consistent with previous ensemble shoreline reconstructions over the same period at Truc Vert beach [START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF], which produced similar shoreline trends and 1998 envelope of shoreline positions, providing a qualitative validation of the wave ensemble used here.

Using the SF model in the 'future scenario', the GSA identified the inherent wave forcing variability as the dominant driver of the modelled shoreline variance, with the estimated Sh's value of the wave chronology roughly doubling the cumulated effects of the ks + and φ parameters over the simulated period Although the Y09 model does not allow to propagate the uncertainties on modelled shoreline towards the past, for the sake of completeness of this analysis, we repeated the 'past scenario' GSA using the Y09 model instead of SF, where the initial condition was set to the mean shoreline position obtained during the model calibration (2012-2019) (Text S3 of Supplementary Material). Unlike the SF 'past scenario' application, the use of Y09 results in no long-term trend of the shoreline ensemble, with seasonal fluctuations around a stable position equal to the mean shoreline position observed during the calibration (2012-2019) (Figure S4). The GSA results show that, on average, the uncertainties in modelled shoreline trajectories are driven by the stochastic wave chronology for 44% while the remaining 56% of shoreline uncertainty is partitioned among the model free parameters (Figure S5).

Future scenarios

The differences in ranges of shoreline projections observed for the Y09+B and SF+B applications reflect the different representations of the physical mechanisms and assumptions underlying the respective disequilibrium formulations. These behaviours were also observed and analysed in previous studies that described the different model responses to short-and long-term variations in the wave forcing (D'Anna et al., 2021a;[START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF]. The SF equilibrium condition is defined by past wave conditions over a fixed time proportional to the beach memory (φ), so that the equilibrium state of the beach adapts dynamically to the wave climate and its temporal variations. Instead, the Y09 equilibrium condition is defined by the previous shoreline position and the fixed equilibrium energy function (Eeq=a1Y+a2), which does not respond to past wave forcing variations beyond a time scale implicitly determined by the fixed value of the a1 and a2 parameters. The different responses to multiannual variations in wave conditions are clearly illustrated by the individual model realizations reported in Figure 5b-e (yellow and green lines). First, the variance of modelled shoreline trajectories is constrained within a limited bound when using the Y09 model, while this is not necessarily the case for SF, consistently with the analysis of [START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF]. However, while characterized by different chronologies, the members of the wave ensemble used here preserve consistent seasonal to multiannual statistics, which somehow constrain the SF envelope of shoreline ensemble to a fixed possible range. The Y09+B and SF+B ensemble shoreline projections also show notable differences in the trends and the ranges of modelled shoreline positions over the initial projection period. In fact, SF equilibrium state over the initial projection period is influenced by the deterministic wave conditions prior to 2021 and the current (random) wave conditions, until the deterministic conditions are 'forgotten'. Given that the two years prior to 2021 are characterized by high-energy winters, the equilibrium condition over the initial projection period is more likely to result in a constructive (accretive) disequilibrium. Over this same period, the Y09 equilibrium condition is only affected by the current (random) wave energy, resulting in an immediate spread of model results and possible trajectories. The narrowing of ensemble ranges observed over the last years of the SF+B projections derives from the assumption that the SF parameter r (Section 2.3.1/Equation 4), which defines the ratio between ks + and ks -, is such that the trend in modelled shoreline is conditioned to the trend of the wave forcing (Splinter et al., 2014b). Therefore, given that all the wave time series are constrained to the same long-term trend, SF will tend to reproduce the same trend by the end of the simulated period for any selected wave time series. The latter model limitation suggests that the assumption of a fixed linear relationship between erosion/accretion response rates may not valid on such long time scales where, instead, the values of model free parameters are likely to be non-stationary [START_REF] Ibaceta | Enhanced Coastal Shoreline Modeling Using an Ensemble Kalman Filter to Include Nonstationarity in Future Wave Climates[END_REF]. The SF assumption on the r parameter has also implications on the GSA, as it drives the SLR's and waves' Sh's to diverge rapidly towards high and low values, respectively (Figure 7d,e), once again, questioning the validity of r's formulation in the context of long-term simulations.

The GSAs for the 'future scenarios' suggest that the uncertainties on wave chronology steadily dominate the variance of shoreline projections until the effects of SLR start concurring to shoreline uncertainties.

The estimated Sh's show that for the simulated scenarios the uncertain future SLR always becomes the primary driver of modelled shoreline variance after 2060. Compared to previous GSA applications on the same study site, but based on deterministic wave time series (D'Anna et al., 2021a), the inclusion of ensemble wave forcing in future shoreline projections postponed the onset of the SLR dominance on the uncertainties of the model results by ~10 years.

Waves uncertainties, interactions and correlation effects

Given the connection between model free parameters and wave forcing in the Y09 and SF models formulations (Splinter et al., 2014b;[START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF], the interaction of these variables and the respective uncertainties within the model may drive higher-order effects contributing to the results uncertainties. Moreover, the correlations among the model parameters of both shoreline models (Figure 4a,b) may further influence the decomposition of the results variance. While a full decomposition of the model results variance that allow discriminating main and interaction terms, and correlation effects is still a matter of ongoing research [START_REF] Mara | Polynomial chaos expansion for sensitivity analysis of model output with dependent inputs[END_REF], here we investigate the presence of interaction and correlation effects by comparing the Sh's with the corresponding individual main effects estimated with first-order Sobol' indices (Si's) (Equation 9, Section 2.4.2).

The difference between Shi and S'i associated to the variable Xi provides an indication of the presence of the cumulated effects of the cumulated terms of correlations and interactions between Xi with the other uncertain inputs. The comparison between Sh's and Si's for the 'past scenario' (Figure 9) suggests the presence of such correlation and interaction effects, between the free model parameters and the wave conditions, on modelled shoreline uncertainties when using the SF model. Indeed, for this scenario the cumulated differences between Sh's and Si's of all the uncertain variables make up ~15% of the total results variance on average, and up to ~75% (near 2012). The latter highlights the risk of underestimating the influence of the model inputs uncertainties on the results uncertainties, and suggests the potential occurrence of contributing waves-parameters interactions. Interestingly, we observe that the seasonal signal of the model sensitivity to the free parameters (higher sensitivity during low-energy periods) is associated to the first-order effects (Figure 9), which is in line with previous assessments of main effects on the modelled shoreline [START_REF] D'anna | Impact of model free parameters and sea-level rise uncertainties on 20-years shoreline hindcast: the case of Truc Vert beach (SW France)[END_REF].

We acknowledge however that the afore-described comparison analysis should only be considered a qualitative indication of the joint impact of interactions and correlations, which might be larger than the differences observed here. This is related to the effects of correlation among uncertain inputs that are partitioned among the correlated variables in the Shapley effects (by construction), but may be accounted redundantly in the Sobol' indices of the correlated variables (Section 2.4.2), causing a possible overestimation of first-order Sobol' indices, hence introducing some biases in the differences Shi-Si.

Dynamic model free parameters

The analysis in Section 5.1 highlighted some shortcomings of both the SF and Y09 approaches: the Y09's limits in reproducing large interannual shoreline patterns, and the SF long-term trend imposed by the fixed r parameter. The respective limitations can be attributed, at least partially, to the assumption of stationary model parameters. Indeed, the SF's free parameters values have been observed to be dependent on time variations of the wave forcing [START_REF] Ibaceta | Enhanced Coastal Shoreline Modeling Using an Ensemble Kalman Filter to Include Nonstationarity in Future Wave Climates[END_REF][START_REF] Splinter | Rapid adjustment of shoreline behavior to changing seasonality of storms: observations and modelling at an open-coast beach[END_REF], while Y09's parameters prevent the equilibrium condition from adapting to shifts in the wave forcing (Vitousek et al., 2021), as discussed in Section 5.1. Therefore, we make a further step to explore the potential impact of the assumption of stationary model free parameters on long-term shoreline projections and run two more sets of 10000 simulations using dynamic values for some of the Y09 and SF parameters for the 'future scenario' RCP8.5. [START_REF] Vitousek | The Application of Ensemble Wave Forcing to Quantify Uncertainty of Shoreline Change Predictions[END_REF] proposed a physical interpretation of the Y09 free parameters, where the a2 parameter in Equation 7 is identified as fixed background wave energy (a2 =𝐸 ̂) that contributes to define the site-specific equilibrium time scale and range of shoreline excursion. Based on site-specific calibration of the model, the optimized value of √ 𝑎 2 was found to approximate that of the average significant wave height (𝐻 𝑠 ̅̅̅ ) over the observed period, so that Equation 7 becomes:

𝐸 𝑒𝑞 (𝑌) = 𝑎 1 𝑌 + 𝑎 2 ≈ 𝑎 1 𝑌 + 𝐻 𝑠 2 ̅̅̅̅ (13)
We run an ensemble of simulations using a modified and calibrated Y09 model (hereon Y09*) based on Equation 13, and we define the background wave energy at a given time as the average significant wave height over the past n years (𝐻 𝑠 ̅̅̅ = 𝐻 𝑠 ̅̅̅ (𝑡)). This assumption not only provides a dynamic character to Y09's equilibrium condition, but also brings Y09 conceptually closer to SF in that the time-window 'n' is analogue to SF's 'beach memory' φ. For this application, we set n=5 years, which correspond roughly to the upper limit of SF's φ values (Table 1). The Y09* model achieves slightly improved skill than Y09 (Section 3.1), with RMSE=4.95 m and R 2 =0.82 on over the calibration. Over the validation period the Y09* model produces an ensemble mean (minimum) RMSE of 8.1 m (3.2 m) and a mean (max) R 2 of 0.89 (0.92), showing a net improvement compared to the Y09 ensemble (Section 4.1). The individual shoreline trajectories (yellow and green lines in Figure 10a) show how Y09* produces interannual fluctuations that were not visible in the Y09 results (Figure 5e). The use of a time-varying a2 parameter produces a more irregular envelope of simulated shoreline positions compared to the results obtained with stationary parameters, with larger amplitudes of the shoreline positions envelope (Figure 10a). It is to be noted that since a2 is fully defined as a function of Hs and the n value is fixed, a2 is no longer an uncertain model parameter and Y09* has less uncertain inputs than Y09.

In the SF model, the ratio of erosion/accretion shoreline response rate (r) is calculated as a function of the trend in wave forcing over the full extent of the simulated period (Equation 4). Here, we relate r to the forcing wave climate dynamically and generate a time series of r values by applying Equation 4 over a progressive time window instead of the full simulation period. At a given time step (t), r (t) is calculated based on the previous n years (n < simulation period) of wave forcing, as follows: The latter applications show that introducing non-stationarity in the Y09 and SF parameters can significantly affect the modelled shoreline response to the wave climate variability and potentially improve the models performances. This encourages future research efforts towards unravelling the inherent relationship between the temporal variability of free parameters and wave climate, and move towards more comprehensive shoreline models. However, these simulations are based on some arbitrary assumptions and are not to be interpreted as an improved assessment of shoreline projections.

𝑟(𝑡) = | ∑ 〈𝐹 + 〉 𝑡 𝑡-𝑛 ∑ 〈𝐹 -〉 𝑡 𝑡-𝑛 | (

Assumptions and limitations

One of the main assumptions of this work is that the current wave climate remains unchanged in the future. In fact, the wave ensemble used here assumes that all the simulated wave time series maintain the same seasonal to multiannual statistics observed in the past wave data, and accounts only for the intrinsic uncertainty associated with the short-term variability of wave conditions. However, climate change could alter the future wave climate in some regions of the world [START_REF] Morim | A global ensemble of ocean wave climate projections from CMIP5-driven models[END_REF], potentially affecting the evolution of sandy shorelines. On one hand, the wave climate over the northeast Atlantic region is not expected to change significantly over the 21 st century [START_REF] Bricheno | Future Wave Conditions of Europe, in Response to High-End Climate Change Scenarios[END_REF][START_REF] Morim | Robustness and uncertainties in global multivariate wind-wave climate projections[END_REF], with dynamically downscaled projections estimating a ~2 mm/year decrease of winter mean Hs in the Truc Vert area [START_REF] Charles | Climate change impact on waves in the Bay of Biscay, France[END_REF]D'Anna et al., 2021a). On the other hand, the limited extent of the hindcast wave data used to generate the wave ensemble (26 years) hinders the inclusion of uncertainties on wave forcing variability on time scales longer than a few decades. In addition, changes in the temporal distribution of extreme events may affect the long-term wave climate cycles [START_REF] Lobeto | Future behavior of wind wave extremes due to climate change[END_REF], and in turn the multi-annual shoreline variability. The inclusion of uncertain future waves in shoreline projections requires wave ensembles accounting for different realistic climate-driven wave climate patterns, such as multi-model ensembles or emulated ensembles based on future climate. To the Author's knowledge, such ensembles do not exist.

While the ensemble of wave time series generated with the wave emulator capture the overall characteristics of the reference past 26-year wave series, the comparison of seasonal wave energy distribution shows some differences (Section 2.4.1). Indeed, the emulator slightly underestimates extreme Hs in the winter, and overestimates them in the remaining seasons. These differences are likely due to a smoothing effect of the seasonality produced by the emulator, together with the strong seasonal variability of the wave climate observed at Truc Vert. However, these differences tend to compensate throughout the different seasons, preserving the overall yearly properties of the wave climate (Figure 2d), so that the simulations are only expected to underestimate seasonal excursions of the modelled shoreline.

In this work, we estimated the future local relative MSL based on SROCC projections. Recently, the IPCC 6 th Assessment Report (AR6) provided updated SLR projections to 2150 [START_REF] Magnan | Cross-Chapter Box 9: Integrative Cross-Chapter Box on Low-Lying Islands and Coasts[END_REF], including new considerations on the Antarctic Ice-sheet dynamic contribution to SLR, and based on combined Shared Socio-economic Pathways (SSP) [START_REF] Riahi | The Shared Socioeconomic Pathways and their energy, land use, and greenhouse gas emissions implications: An overview[END_REF] and RCP scenarios. However, our sea level projections (2021-2100) remain close to the medium confidence projections published in AR6 Working Group I -Chapter 9 (within the order of 1 cm difference).

Finally, we modelled SLR-driven shoreline recession using the Bruun model, which builds on a number of assumptions that limit its applicability to a reduced selection of beaches [START_REF] Cooper | Sandy beaches can survive sea-level rise[END_REF]. As

Truc Vert beach is a cross-shore transport dominated, uninterrupted beach with large accommodation space most of the Bruun model assumptions are satisfied. However, alternative rule-based approaches [START_REF] Mccarroll | A rules-based shoreface translation and sediment budgeting tool for estimating coastal change: ShoreTrans[END_REF], or dynamically coupled wave-driven and sea-level equilibrium approaches (D'Anna et al., 2021b) could be implemented for applications of this framework to a wider range of sites. In addition, here the Bruun Rule combined linearly with Y09 and SF with no feedbacks between wave-and SLR-driven modelled processes, preventing any interaction between wave forcing and SLR uncertainties. Approaches that combine wave-and sea-level-driven processes in a coupled fashion (i.e.

with feedback between the models), may reveal further contributions of the wave forcing uncertainties to modelled shoreline change.

Conclusions

We applied a Global Sensitivity Analysis to long-term shoreline evolution modelling at the cross-shore transport dominated Truc Vert beach (France) including the intrinsic uncertainties associated with unknown wave chronology and using two different equilibrium shoreline models. Our results suggest that uncertainties in wave chronology play an exceedingly important role in probabilistic shoreline modelling. When modelling shoreline change with no uncertainties on SLR (e.g. 'past scenario' and near-future projections), the uncertain wave chronology is responsible for up to 83% of the modelled shoreline variance. We also find that the interactions between the uncertain wave chronology and model parameters within the models and the correlation among the model parameters can significantly influence shoreline predictions uncertainties, highlighting the interplay between intrinsic and epistemic uncertainties. The future shoreline projections suggest that the onset of the SLR dominance on the uncertainties of the model results can be expected to start from 2060, with strong implications for future research priorities and decision-making in climate change adaptation planning where coastal settings are similar to Truc Vert beach. While the two equilibrium modelling approaches compared in our applications (ShoreFor and Yates) show different responses to wave climate variability and its uncertainties, the main conclusions drawn above are valid for both approaches. We also found that applying non-stationary model free parameters by linking them to the wave climate variability can significantly alter the model behaviour, and further research efforts in this direction could improve the understanding of model uncertainties.

Figure 1 .

 1 Figure 1. (a) Map of the southwest coast of France with the location of Truc Vert beach (green), and wave hindcast grid point co-located with the CANDHIS in situ wave buoy (red); Photographs of Truc Vert beach and dune landscape from (b) UAV aerial view (photo by B. Castelle), and (d) foredune view (photo by S. Bujan); (c) 4 km alongshore-averaged beach-dune profile from merged 2008 topo-bathymetry (submerged beach) and 2018 UAVphotogrammetry digital elevation model (emerged beach and dune); (e) Time series of alongshore-averaged shoreline (1.5-m elevation proxy) position from April-2005 to June-2021 derived from the topographic surveys. Vertical dashed lines divide the data prior to March-2012 (cian circles) from data used for model calibration (2012-2019) with the respective linear trend and validation (2019-2021) (white circles) (see Section 3.1).

  model (herein Y09) the disequilibrium condition is expressed in terms of incident wave energy (ΔE = Eeq -E), where the equilibrium state Eeq is defined as a function of the current shoreline position Y. In Y09, the rate of shoreline change is calculated as follows: 𝑑𝑌 𝑑𝑡 = 𝑘 𝑦 +/-𝐸 0.5 (𝐸 𝑒𝑞 (𝑌) -𝐸 ) (6) where ky +/-(m² s -1 /m) is the response rate parameter, E(m²) is the wave energy, Y(m) is the current shoreline position, and Eeq (Y) is the wave energy that would cause no change to the current shoreline position Y. Similarly to SF, the response rate parameter ky +/-indicates the rate (or efficiency) of shoreline response to the wave forcing, and associates different values for accretion (ky + , ΔE > 0) and erosion (ky -

Figure 2 .

 2 Figure 2. (a) Time series of Hs (black line) and 3-month averaged H²sTp (red line) from 1994 to 2021 from NORGAS-UG hindcast model; Envelope (min. and max.) (red area) and mean ± σ interval (yellow area) of the 200 (b) Hs and (c) 3-month averaged Hs 2 Tp time series from 2021 to 2100, generated using the wave emulator, with mean Hs (orange line) and the 3-month averaged Hs 2 Tp of an individual time series; quantile-quantile comparison between the Hs of hindcast data and each member of the wave ensemble for (d) the full time series, and (e) their seasonal repartition.

Figure 3 .

 3 Figure 3. Method outline of a generalized case (black box) and the Truc Vert beach applications (red box) for past and future different sea-level rise and shoreline modelling scenarios. Figure modified from D'Anna et al. (2020).

Figure 4 .

 4 Figure 4. Joined distributions of (a) the ShoreFor (ks + , φ) and (b) Yates (ky +/-, a1, a2) free parameters with the corresponding Pearson correlation coefficients ρ; and probability distributions of 2021-2100 sea-level rise projections for (c) the RCP4.5 and (d) RCP8.5 scenarios, with the respective 66% confidence interval (dark-shaded areas) and 5 th -95 th percentile range (light-shaded area).

  Following D'Anna et al. (2021a), we produced probabilistic MSL projections assuming a Gaussian distribution (in line with[START_REF] Hunter | Towards a global regionally varying allowance for sea-level rise[END_REF], with median and standard deviation corresponding to the yearly median and likely range SROCC estimates for the RCP4.5 and RCP8.5 scenarios (Figure4c,d). As the SROCC projections refer to the year 2007, we re-fitted the projections between 2021 and 2100, conditioning the MSL median and likely range to the corresponding values observed in 2021. See Text S1 and FigureS3of Supplementary Material for more details on this correction. As the characteristics of the future MSL probability distribution are timedependent, in the following applications MSL time series were defined for a given percentile (p) by extracting the sea level corresponding to p at each year.

  Figure 5a-c and Figure 5d,e show the envelopes (min-max) and likely ranges (µ±σ) of shoreline trajectories obtained using the SF and Y09 approaches, respectively, for the different scenarios. All modelled ensembles consistently fit the shoreline data over the calibration period (2012-2019) and well reproduce the large erosion associated with the particularly energetic winter of 2013-2014. The models also reproduce the shoreline variations observed during the validation period (2019-2021), with a mean determination coefficient R 2 =0.73 and a mean RMSE=10.65 m. Over the validation period, the modelled shoreline ensembles capture the erosive trend associated with the increasing interannual winter wave energy in 2019-2021 (Figure

- 7 .

 7 6 m] ([-59.2 m; 46.9 m]) in 1998 to[-8.3 m; -1 m] ([-15.1 m; 6.3 m]) in 2012. The application of ensemble wave forcing(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) results in a 'richer' envelope of possible shoreline trajectories and visibly larger prediction variance compared to the calibration and validation period (2012-2021) where a deterministic wave time series was applied (Figure5a). The ensemble results also envelope the shoreline data prior to 2012, which was discarded from the model calibration but provides a realistic reference of shoreline variability.

Figure 5b -

 5b Figure5b-e show the 'future scenarios' of ensemble shoreline projections modelled combining the Bruun model with SF and Y09, respectively, for two SLR scenarios (RCP4.5 and 8.5). All models were forced with the same random selection of wave time series, and no feedbacks were allowed between the wave-and SLR-driven models. Therefore, any long-term shoreline trend induced by each wave series of the ensemble is reproduced identically for the RCP4.5 and 8.5. Hence, the results obtained for the different RCP scenarios only differ in the magnitudes of the SLR-driven erosional trends. Indeed, using either SF or Y09, the RCP8.5 scenario produces a more eroded shoreline position in 2100 and a larger variance compared to the RCP4.5 scenario. Likely ranges and envelope values of the 2100 simulated shoreline position are summarized in Table3for each 'future scenario', including the corresponding ranges resulting from the application of the Bruun Rule alone.

Figure 5 .

 5 Figure 5. Ensemble shoreline evolution simulated using the (a-c) ShoreFor and (d,e) Yates, combined with Bruun models, including the envelope (light-shaded areas) and the standard deviation intervals (µ±σ) (dark-shaded areas) of the 10000 simulated shoreline trajectories, respectively, over (a) the 1998-2021 period ('past scenario'), and from 2012 to 2100 ('future scenarios') for the (b,d) RCP4.5 and (c,e) RCP8.5 future sea-level rise scenarios. Yellow and green lines (b-e) are individual model realizations within the ensemble. Black dashed lines delimit the calibration and validation periods, including shoreline observations (white circles). Cian-coloured circles (a) indicate shoreline data prior to March 2012.Table3. Model results at 2090 and 2100 obtained for the four 'future scenarios' and the applications of the Bruun Rule alone, including the likely range (µ±σ) and envelope (min-max) of simulated shoreline position, and the mean of 1-year averaged shoreline trajectories.

Figure 6 .

 6 Figure 6. Global Sensitivity Analysis results for the 'past scenario', including time series of: (a) the ensemble simulated shoreline positions from 1998 to 2012; and Shapley Effects of the ShoreFor parameters (b) ks + and (c) φ, and (d) wave chronology, with the respective linear fit (black lines).

Figure 7 .

 7 Figure 7. Global Sensitivity Analysis results for 'future scenarios' simulated using ShoreFor for the RCP4.5 (blue) and RCP8.5 (orange) SLR scenario. Time series of (a) standard deviation of simulated 2021-2100 shoreline positions, and Shapley Effects of model parameters (b) ks + and (c) φ, (d) sea-level rise rate, and (e) wave chronology.

Figure 8 .

 8 Figure 8. Global Sensitivity Analysis results for 'future scenarios' simulated using Yates' model for the RCP4.5 (blue) and RCP8.5 (orange) SLR scenario. Time series of (a) standard deviation of simulated 2021-2100 shoreline positions, and Shapley Effects of model parameters (b) a1 (c) a2, (d) ky + , (e) ky -, (f) sea-level rise rate, and (g) wave chronology.

(

  Figure 7b,c,e). All Sh's time series exhibit a very mild seasonal signal between 2005 and 2012, associated with higher model sensitivity to the ks + and φ parameters during low-energy periods (summers), and stabilize over time as the model variance increases (Figure 7b,c). Compared to previous GSAs where only model free parameters governed the modelled shoreline variance (D'Anna et al., 2020), the effects of uncertainties in wave forcing variability override and buffer the complex seasonal patterns of the model sensitivity to model free parameters. Our findings are also in line with the analytical results of Vitousek et al. (2021), who showed that the intrinsic uncertainties on wave chronology can more than double the contribution of the epistemic uncertainty on calibrated model parameters to the shoreline ensemble variance.

Figure 9 .

 9 Figure 9. Comparison of Shapley Effects (green lines) and 1 st order Sobol' (black lines) sensitivity indices (s.i.) of the ShoreFor free parameters (a) ks + , (b) φ, and (c) wave chronology, for the 'past scenario' ensemble.

  14)Such formulation preserves the original SF concept that links the trends of shoreline change and wave climate but allows different multi-decadal wave climate patterns to result in different shoreline positions by the end of the simulated period. We calibrated a modified SF (hereon SF*) model based on Equation14with n= 5 years and used it to run a new ensemble of simulations. The SF* model produces an RMSE=6.8 m and a R 2 =0.71 on the simulated period, showing some improvement compared to SF (Section 3.1). The model skill also improved over the validation period, with a mean (min) RMSE of 8.62 m (3.05 m) and a mean (max) R 2 of 0.85 (0.93). The use of non-stationary values for the r parameter results in a progressive increase of the modelled shoreline variance, due to the increasing number of possible sequences of high-and low-energy winters (Figure10b). The values of erosion/accretion rate r computed with Equation 14 occasionally attains values above 10 and near 0. This is likely to cause the σ (envelope) ranges of the modelled shoreline to spread over time up to roughly [-10 m; -110 m] ([150 m; -250 m]). However, while the SF* tested here is based on some arbitrary assumptions (i.e. n=5, no limits in the range of r variability), the model produces different long-term trends in response to the different wave forcing time series, and do not converge to the same final shoreline position at the end of the simulated period.

Figure 10 .

 10 Figure 10. Ensemble of simulated shoreline trajectories from 2021 to 2100 for the RCP8.5 future sea-level rise scenario using the Bruun Rule with the modified (a) Yates et al. (2009) and (b) ShoreFor models, including: the envelope (light-shaded areas) and the standard deviation intervals (µ±σ) (dark-shaded areas); two individual model realizations (yellow and green lines). Black solid lines indicate the limits of the envelope of the corresponding simulations obtained using the original Yates et al. (2009) and ShoreFor models.

Table 1 .

 1 Optimised combinations of cross-shore model free parameters, and respective range of variation in the probability distributions.

	Model	Model parameter Optimised value Distribution range
	ShoreFor (SF)	ks + [m 1.5 s -1 W -0.5 ]	5.0 x10 -8	[1.6; 9.8] x10 -8

Table 2 .

 2 Past and future simulated ensembles based on 10000 different combinations of model parameters, SLR percentile and wave time series, for two future Representative Concentration Pathways (RCP) and two different wave-driven modelling approaches with the Bruun Rule.

	Scenario SLR scenario	SLR-driven shoreline model	Wave-driven shoreline model	Dataset name	Combinations	Simulated period
	Past	Reconstructed	Bruun Rule	ShoreFor (SF)	SF+B/Past	10000	1998-2021
	Future	RCP 4.5 RCP 8.5	Bruun Rule Bruun Rule	ShoreFor (SF) Yates (Y09) ShoreFor (SF) Yates (Y09)	SF+B/RCP4.5 Y09+B/RCP4.5 SF+B/RCP8.5 Y09+B/RCP8.5	10000 10000 10000 10000	2012-2100

Acknowledgements

This study includes the monitoring study site of Truc Vert labeled by the Service National d'Observation (SNO) Dynalit (https://www.dynalit.fr) with additional support from Observatoire Aquitaine de l'Univers (OASU) and Observatoire de la Côte Aquitaine (OCA). The authors thank G. Le Cozannet and R. Thiéblemont involved in the sea level data and early discussions; S. Bujan, S. Ferreira and V.

Marieu involved in the topographic data; SONEL for past sea levels and vertical land motion data (https://www.sonel.org/-GPS-.html); the Integrated Climate Data Center at the University of Hamburg (https://icdc.cen.uni-hamburg.de/en/ar5-slr.html) for SLR projections data; and the developers of the R-Studio sensitivity package used here to evaluate the Shapley Effects.

Funding: This work is co-financed by the Make Our Planet Great Again (MOPGA) national program Grant [927923G] and the BRGM (French Geological Survey). BC is funded by Agence Nationale de la Recherche (ANR) Grant [ANR-17-CE01-0014]. LC and FJM acknowledge the funding from the Spanish Ministry of Science and Innovation, project Beach4cast [PID2019-107053RB-I00].

Data Availability Statement

All data used and produced within the current work is available at the following Mendeley data online repository: DOI: 10.17632/cnwyhh475s.1