Which data assimilation method to use and when: unlocking the potential of observations in shoreline modelling

Shoreline predictions are essential for coastal management. In this era of increasing amounts of data from different sources, it is imperative to use observations to ensure the reliability of shoreline forecasts. Data assimilation has emerged as a powerful tool to bridge the gap between episodic and imprecise spatiotemporal observations and the incomplete mathematical equations describing the physics of coastal dynamics. This research seeks to maximize this potential by assessing the effectiveness of different data assimilation algorithms considering different observational data characteristics and initial system knowledge to guide shoreline models towards delivering results as close as possible to the real world. Two statistical algorithms (stochastic ensemble and extended Kalman filters) and one variational algorithm (4D-Var) are incorporated into an equilibrium cross-shore model and a one-line longshore model. A twin experimental procedure is conducted to determine the observation requirements for these assimilation algorithms in terms of accuracy, length of the data collection campaign and sampling frequency. Similarly, the initial system knowledge needed and the ability of the assimilation methods to track the system nonstationarity are evaluated under synthetic scenarios. The results indicate that with noisy observations, the Kalman filter variants outperform 4D-Var. However, 4D-Var is less restrictive in terms of initial system knowledge and tracks nonstationary parametrizations more accurately for cross-shore processes. The findings are demonstrated at two real beaches governed by different processes with different data sources used for calibration. In this contribution, the coastal processes assimilated thus far in shoreline modelling are extended, the 4D-Var algorithm is applied for the first time in the field of shoreline modelling, and guidelines on which assimilation method can be most beneficial in terms of the available observational data and system knowledge are provided.


Introduction
Sandy shorelines are among the most dynamic coastal environments on Earth, cover one-third of icefree coasts, host key socioeconomic activities and provide valuable flood protection for many coastal cities worldwide (Toimil et al 2023).Predicting how they react to natural and anthropogenic drivers is becoming one of the most important challenges in coastal science (Toimil et al 2020a, Splinter andCoco 2021).Physics-based modelling is the most viable alternative for efficiently generating various related predictions, from short-term cross-shore sediment transport (Miller and Dean 2004, Yates et al 2009, Splinter et al 2014) to long-term shoreline changes driven by longshore transport (Pelnard-Considère 1956) and water level changes (Bruun 1962, McCarroll et al 2021).
For shoreline models to provide reasonable predictions, their free parameters need to be tuned via calibration (Splinter et al 2013, Montaño et al 2020), which presents several challenges.One challenge is the computational cost associated with traditional calibration techniques, including Monte Carlo variants such as general likelihood uncertainty estimation (Simmons et al 2017) or minimization of a cost function using gradient-free algorithms (Hewageegana and Canestrelli 2021).
Another challenge arises from the fact the validation of physics-based shoreline models is limited to a few monitored coastal sites (Turner et al 2016, Ludka et al 2019, Castelle et al 2020).In this sense, the release of satellite imagery and its use in coastal science (Luijendijk et al 2018, Vos et al 2019, Sabour et al 2020, Sánchez-García et al 2020), as well as the emergence of new forms of monitoring such as citizen initiatives (Harley and Kinsela 2022) and the application of the Internet of Things (IoT) technology (Tien et al 2023), have enabled important steps forwards.However, guidelines for determining the quantity (the length of the observation campaign and the sampling frequency) and quality (the accuracy of the measurements) of the data required for calibration to provide reliable estimates are still lacking.
Increased data availability has boosted the use of artificial intelligence (AI) and data assimilation (DA) in geophysical modelling in recent years (Kim et al 2021, Sonnewald et al 2021, Kelp et al 2023).While AI can be used to infer system dynamics directly from observational data (Kuglitsch et al 2023), DA blends a dynamic model with observations to better estimate the true state of a system and eventually tune its unknown parameters (Asch et al 2016).Unlike traditional calibration methods (Simmons et al 2017, Hewageegana andCastrelli 2021), DA allows for the joint tuning of parameters and predictions so that the true (unknown) state of the system can be efficiently reproduced considering the uncertainty of the model and the observations, through just one or a reduced number of model runs.There are two main families of DA methods: statistical and variational methods.The former requires the minimization of the uncertainty of the model error (Kalman 1960), while the latter seeks to minimize a cost function that represents the mismatch between the model and the observations (LeDimet and Talagrand 1986).
Variational methods, such as 4D-Var (the 4dimensional variational method), and statistical methods, such as Kalman filter-based algorithms, have both been widely used for data assimilation in meteorology and oceanography (e.g.Houtekamer et al 2005, Neveu et al 2016, Stopa 2018).The choice between them boils down to a trade-off between accuracy, computational complexity, the characteristics of the observations and the characteristics of the type of analysis to be conducted (Kalnay et al 2007).Kalman filters are more flexible and easier to implement than 4D-Var in systems with real-time data updates (Lorenc 2003).Nevertheless, despite its complexity and the need for an adjoint model, 4D-Var also has advantages.It considers the entire time interval of interest as a whole to optimize the estimation of the initial state.This can be very useful in systems where the observations are irregularly distributed over time.Additionally, 4D-Var can be more accurate for nonlinear and dynamic systems (Caya et al 2005).In the field of coastal impacts, the latter advantage can translate into improvements in addressing changing trends and nonstationarities such as those arising from climate change.
To date, there have been few DA applications in shoreline modelling, and these applications have been limited to statistical methods (Kalman filters) combined with observations, whereas variational algorithms have remained unexplored.Shoreline prediction with DA was first introduced by Long and Plant (2012) by combining the extended Kalman filter (eKf) with a cross-shore model plus a linear term representing longshore processes using synthetic observations.Vitousek et al (2017) developed an integrated longshore and cross-shore model in combination with an eKf fed with LiDAR observations, while Alvarez-Cuesta et al (2021) used satellite data together with an eKf and a multiprocess shoreline model.Whereas the LiDAR observations were accurate (with a root mean square error, RMSE, of <1 m) and recorded once a year, the satellite observations were noisy (RMSE ∼ 10 m) but sampled more frequently (biweekly).Ibaceta et al (2020) employed a crossshore model and an Ensemble Kalman filter (EnKf) to track nonstationary parameters using observations from video cameras (RMSE ∼ 0.5 m, sampled on a daily to weekly basis).Vitousek et al (2021) used an EnKf to assess the uncertainty of shoreline predictions due to wave forcing while additionally using cameraderived observations.Recently, Vitousek et al (2023) introduced a localized EnKf together with a previously developed multiprocess shoreline model that accounted for spatial dependencies in the littoral cells of California's coasts.They concluded that the use of satellite-derived shorelines yielded results similar to those obtained using monthly in situ GNSS data (RMSE ∼ 1 m) from Ocean Beach.
The free model parameters to be tuned through DA should not change with variations in climate forcing (waves and water levels) because, by definition, they are independent of any other model input.In reality, however, this expectation may not hold (Ibaceta et al 2020(Ibaceta et al , 2022)), and climate variability may consequently compromise the quality of the calibration process, due to hidden dependencies between the wave climate and the free model parameters that should be disentangled by the assimilation algorithm.Thus, shoreline forecasts may benefit from the incorporation of climate-dependent parametrizations, as the world's coasts are expected to be exposed to unprecedented drivers due to climate change (Toimil et al 2020b, Lobeto et al 2021).
The main objective of this paper is to provide general guidelines for exploiting current shoreline observations to guide physics-based models of shoreline evolution at different scales using the most appropriate assimilation algorithm.To achieve this goal, the coastal processes considered in data-assimilated shoreline modelling are extended; the 4D-Var assimilation algorithm is applied, for the first time to the authors' knowledge, to shoreline modelling; and the data and system knowledge requirements for the success of statistical and variational algorithms are analysed.Conclusions are drawn on the basis of synthetic exploratory cases and are subsequently verified at two real beaches subject to different coastal processes and with different data used for calibration.

Numerical models
In this work, shoreline evolution is regarded as a linear sum of cross-shore and longshore processes.The short-term coastal response to changes in wave action and water levels results in cross-shore sediment transport that moves sediment from the beach face to the offshore bar (erosion) and vice versa (accretion).Conversely, the obliqueness of the wave angle with respect to the coastline drives longshore sediment transport, whose gradients generate shoreline changes on a yearly time scale.
Using a physics-based modelling approach, crossshore processes are modelled following Miller and Dean (2004), as shown in equation ( 1), but in addition, the astronomical tide is included, as in Toimil et al (2017).The model translates changes in waves and water levels to shoreline changes driven by the imbalance between the equilibrium shoreline Y eq c and the instantaneous shoreline Y c where is the accretion/erosion cross-shore rate and Y eq c may be expressed as: ∆y 0 is an empirical parameter that ensures that short-term fluctuations oscillate around a baseline position, W * b is the breaking line width, H sb is the significant breaking wave height, SS is the storm surge, AT is the astronomical tide amplitude, B is the beach profile berm height and h b is the breaking depth.
Longshore changes are modelled at fixed shorenormal transects, which define the cell boundaries where the longshore transport gradients are evaluated.Gradients in longshore sediment transport generate landward or seaward displacements of the active part of the beach profile: Y l is the distance from the onshore point of the transect to the shoreline, d c is the closure depth, ∂x is the transect spacing and Q is the longshore sediment transport calculated according to the Coastal Engineering Research Center (CERC) expression: where K l is an empirical constant and α is the angle between the wave crests (θ b ) and the shoreline (θ s ).

Data assimilation algorithms
Given a set of observations y that are scattered in time and space, DA attempts to find the optimal state x defined by the variables evolved by the dynamic model.Shoreline observations y can be acquired through various means, such as video imagery, satellite imagery or in situ surveys conducted at sampling intervals spanning from hours to several months.These digitalized shorelines, whose spatial resolution depends on the acquisition method, are then interpolated onto a set of model transects.The state x is formed by the model predictions (Y c and Y l for crossshore and longshore processes, respectively) and augmented by the free model parameters (K +/− and K l for cross-shore and longshore processes, respectively).For cross-shore processes, a single transect is considered in the state vector.In contrast, the state vector for longshore transport is formed by Y l and K l at every transect used to discretize the shoreline.
Considering the state and the observations as random variables, the DA problem consists of solving the Bayes equation that gives the posterior distribution of the state (x) given the observations based on the prior state distribution and the likelihood of the observations: where p (y) is the marginal probability density function of the observations.It acts as a normalizing constant that ensures that p (x|y) is a valid probability distribution (integrates to one).A common practise in statistical and variational methods is to assume that the distributions are normal so that they can be modelled in terms of the mean and the covariance matrix (B): , where the superscripts b and a stand for background and analysis, respectively.The initialization of B b is key for the success of DA algorithms (Bannister 2008a); interested readers can refer to the supplementary material for further details.R is the observation error covariance matrix.It is common to assume unbiased and uncorrelated errors (i.e. a diagonal R matrix), where each term is the variance of the observation error distribution function, which depends on the accuracy of the data source (e.g.satellite images, LiDAR, GNSS, video camera footage), plus the variance of the interpolation errors that arise when mapping the temporally and spatially scattered observations to the modelled state (negligible in shoreline modelling).
The covariance matrix for a Kalman filter is denoted by P, and its initial value coincides with the initial background covariance matrix, P 0 = B b 0 .The initial background covariance matrix is constructed by assuming no correlations between the model state and the parameters (Evensen et al 2022).The background covariance matrix is formed from the errors of the state with respect to the background.The error in the initial shoreline position is assumed to be similar to the error in the observations, and the errors on the parameters are of the same order of magnitude as but smaller than the values of the initial constants.The interested reader is referred to the supplementary material for further details on the initialization of P 0 .
In sequential assimilation, there are two computational steps: forecasting, which involves a free model run, and correction, which involves an information flow from the observations to the model.
During the forecasting step (superscript f), the analysis state x a i is propagated in time considering the model's dynamic equation: The eKf propagates the analysis covariance matrix P a i analytically using the tangent linear model M i of M at instant i considering the covariance of the model's uncertainty or process noise Q: During the correction step, the Kalman equations yield the mean x a i+1 and the covariance P a i+1 of the posterior distribution in a sequential fashion every time a new observation y i+1 is available (figure 1(a)): )) H is a nonlinear observation operator that maps the model state to the observation space, and H is its linearization.It is generally a temporal and spatial interpolation that enables the comparison of the model state (i.e. the shoreline) to the observations.I is the identity matrix, and K is the Kalman gain, which weights the model state with respect to the observations.Due to the size of P f i +1 , which scales with the size of the augmented state, this matrix is unfeasible to calculate in many environmental applications.Thus, it can be approximated using a sample of evolved states, leading to the EnKf (Evensen 1994).
In the SEnKf, each ensemble member is analytically evolved following equation ( 6), and the covariance matrix P f i +1 is directly obtained from the ensemble without analytical derivations (figure 1(b)): where X f i +1 is the perturbation matrix defined by subtracting the ensemble mean from each ensemble member's state and N is the number of ensemble members.
In the SEnKf, the observations need to be perturbed for every ensemble member and updated so that the variance according to the SEnKf analysis correctly represents the variance according to the eKf analysis (i.e.(I − K i+1 H) P f i+1 ).Thus, equation ( 8) becomes: )) , Adding noise to the observations via random sampling from N (0, R i +1 ) avoids the divergence of the SEnKf by adjusting the analysis variance.As the magnitude of the sampling noise depends on the error of the observations via the covariance matrix R i +1 , the larger the data source error is, the larger the sampling error introduced in the SEnKf.This fact limits the application of the SEnKf when the variations in the time of the state are relatively small compared to the data source error.
The magnitude of the sampling noise depends on the error of the observations via the covariance matrix R i +1 .
In variational methods, the sought solution is the one that maximizes the posterior probability of p (x|y) (i.e. the mode), and no uncertainty estimate of the analysis state is obtained.In 4D-Var, DA is performed in a time window that enables future observations to be considered in prediction at a previous time step (figure 1(c)).
Considering that inside the time window, the state trajectory is governed solely by the model dynamics (strong-constraint 4D-Var), the minimization of the cost function (equation ( 11)) yields the most probable state where x 0 b is the a priori background state, B 0 is the initial background error covariance matrix, x 0 is the initial state vector, R i is the observation error covariance matrix, K is the number of steps in the time window and x i is the state at time t i of observation y i .Note that ultimately, J is a function of x 0 , as Additionally, if there are no observations at a given step i, Minimizing J requires the term-by-term calculation of the gradient of the cost function to guide a descent algorithm: Solving for ∇J (x 0 ) implies the backpropagation of the innovation vector y i − H (x i ) using the transpose (M T ) or adjoint of the model operator dynamics.Deriving the adjoint of the forward model enables the computation of the gradient of the cost function, as shown in the supplementary material.
The 4D-Var algorithm consists of iteratively running the forward model (equation ( 6)), calculating the cost function (equation ( 11)) and its gradient (equation ( 12)), and applying a descent algorithm to find the optimal initial state vector x 0 .J (x 0 ) is usually nonquadratic, preventing the descent algorithm from reaching the global minimum, which corresponds to the usually unknown ground-truth (reference) solution.Hence, incremental 4D-Var and control variable transformations are usually required in practical applications.These techniques are comprehensively described in the supplementary material.

Description of the test cases
The parameter estimation skill of the different DA algorithms is tested considering the observation data requirements in terms of accuracy, length and sampling frequency and also the required system knowledge in terms of initial conditions and parameter nonstationarity.Cross-shore and longshore processes are studied in isolation following the twin experimental procedure.This method consists of generating a ground-truth simulation (reference) using a set of parameters by freely running the dynamical model.This simulation is then used to generate the observations for other simulations, in which the initial parameter set is altered but the dynamical model is enhanced by the DA algorithms.
Cross-shore processes are modelled at a single transect.The reference cross-shore shoreline position is Y st,0 = 0, and the accretion and erosion constants are K + = 0.002 and K + = 0.04, respectively.Waves and water levels correspond to real time series obtained from Alvarez-Cuesta et al (2023).The equilibrium profile for the cross-shore model is characterized by Dean's parameter (Dean 1991) with a value of A = 0.25.Longshore processes are modelled by considering a 3 km embayment discretized with 15 shorenormal transects.The reference initial shoreline is a parabola calculated such that the cross-shore length of the embayment is 100 m.Null flow conditions are imposed at the boundaries to simulate a closed system, and the wave dynamics are constant, with H sb = 2 m and θ b = −30 • .A uniform longshore constant K l = 50m 1/2 /day is considered at every transect.Further details about the reference state of the system and the initialization of the perturbed simulations for the twin experiments are provided in the supplementary material.
The length of the simulations is 30 years, divided into a calibration phase (DA-guided run) and a validation phase (free run without DA).When analysing the observations and the initial system knowledge requirements, the validation phase covers the last 10 years.When analysing the nonstationary evolution of the system, because the ability of the DA algorithm to track changes in the model parameters from the observations is measured, DA is never switched off, and the calibration period covers the whole time span.The main characteristics of the terms influencing the physics-based process equations (equations ( 2)-( 4)) in the test cases and the specific characteristics of the different DA methods (i.e. the number of ensemble members for the SEnKf and its initialization, the temporal window for 4D-Var, and the definitions of the state covariance matrix B and the process noise matrix Q) are described in the supplementary material.
Different error metrics are used to test the capabilities of the DA algorithms throughout the analysis during the validation phase: • DA benefit (DAB).This metric measures the improvement of a DA algorithm with respect to the background simulation (without DA) in terms of the root mean squared error (RMSE) between the altered simulation and the reference simulation during the validation period: where the subindex alt stands for 'altered' and b for 'background' .
• Parameter error (eK).This metric measures the difference between the DA value at the end of the calibration period and the reference value The subindex r refers to the reference simulation.For longshore processes, eK is calculated based on the mean DA parameter across all transects.• Computational burden (CPU).This quantity is measured in seconds.• Percentage of time within acceptable limits (PTWL).This metric is used to evaluate whether a DA algorithm is capable of tracking nonstationary variations in the model parameters.The acceptable limits are defined as a 20% range of a parameter around its reference value.

The role of the observations: quantity and quality
The performance of the DA algorithms is tested by using observations with different levels of added noise (N ( 0, ϵ 2 ) ), where ϵ includes the error on the observations, as well as variations due to different observation lengths and different sampling frequencies.In figure 2(a), the performance metrics are analysed for different observation lengths and added noise levels considering a fixed sampling frequency of 15 days.In figure 2(b), the sampling frequency is analysed together with the length of the observations for an observation error of 1 m.
This analysis shows that the Kalman filters perform better than 4D-Var for both cross-shore and longshore processes when the reference parameters do not vary in time.For cross-shore processes, the eKf slightly outperforms the SEnKf, as it achieves a DAB larger than 80% in almost every scenario at a smaller computational expense.In the case of longshore processes, the SEnKf outperforms the other algorithms in terms of DAB and eK when the error on the observations is greater than 7 m and when the length of the observations is shorter than 3 years.Conversely, 4D-Var performs well when the error on the observations is less than 4 m for both cross-shore and longshore processes.For greater observation errors, increasing the length of the observations is not clearly correlated with an improvement in DAB.This fact is explained by the random nature of the observations together with the lack of flow dependency of the state covariance matrix between assimilation windows (interested readers are referred to the supplementary material for a detailed explanation).For longshore processes, the length of the observations should exceed 5 years to yield a relative error on the model constant of less than 20% of the reference value of the parameter.From a computational perspective, the SEnKf is the most demanding algorithm because of the use of an ensemble, and its cost increases with the number of observations for simulations exceeding 50 s.The most efficient algorithm in the case of crossshore processes is the eKf, with an average simulation cost of less than 10 s, while 4D-Var is also computationally efficient but becomes more expensive as the length of the observations increases, as the minimization algorithm needs to be executed more times.
Regarding the sampling interval of the observations, the DA algorithms do not show great differences in the 1-60 days range.The SEnKf yields a DAB greater than 80% for all sampling intervals from 1 to 60 days, and the relative parameter error is always smaller than 20% for both longshore and cross-shore processes.The eKf yields a DAB of more than 80% in every case, except when the sampling interval is 60 days and the length of the observations is one year for cross-shore processes or two years for longshore processes.However, the longer the sampling interval is, the more observations are required to reduce the relative error of the parameters from the 20%-40% range to the 0%-20% range.Similarly, 4D-Var yields a DAB of more than 80% in every case for cross-shore processes, but its parameter convergence is compromised when the sampling interval of the observations is greater than 30 days for campaign durations of less than four years.In the case of longshore processes, to guarantee good convergence of the algorithm, observations over more than 3 years are required if the sampling interval is greater than one month.

The role of system knowledge: prior knowledge and temporal evolution (parameter nonstationarity)
The initial state of the perturbed simulations is modified by taking different quantiles from the upper and lower tails of the CDF of the reference simulation, as detailed in the supplementary material.Observations with no added noise are assimilated every 15 days.
In figure 3(a), the performance metrics are analysed for different campaign durations and initialization values.The Kalman filters require more precise initialization than does 4D-Var for cross-shore processes, which is especially noticeable for extreme quantiles (0.999).In this case, the SEnKf requires more than 10 years of observations; the eKf, 5 years; and 4D-Var, 2 years to achieve a DAB greater than 80%.For longshore processes, the Kalman filters converge to the reference value even for extreme quantiles, but 4D-Var requires more than 5 years to successfully converge in the case of the 0.999 quantile.The different time scales associated with longshore and cross-shore processes may explain the different convergence rates.
The skill of the DA algorithms in tracking the variations in time of the free model parameters (parameter nonstationarity) is measured by forcing the cross-shore and longshore numerical models using square waves.A square wave is formed by pulses such that the value of the varying parameter alternates between two discrete levels in a regular and repetitive manner.The pulse period, defined as the time it takes for one complete cycle of the waveform to occur, is used in combination with the error of the observations to rank the different DA algorithms.The DA performance indicators are calculated for the whole simulation period, as the calibration and validation phases coincide, and biweekly observations are used for assimilation.
In figure 3(b), the performance of the DA methods is analysed for different pulse periods of the square wave and for different observational errors.In this case, 4D-Var clearly outperforms the Kalman filters for cross-shore processes.All DA algorithms perform better for longer pulse periods, but the error on the observations is found to be a critical parameter when analysing parameter nonstationary.For cross-shore processes, the SEnKf never reaches a DAB greater than 80%, even for clean observations.The eKf performs slightly better, but clean observations are required for this algorithm to show an improvement with respect to the background simulation.In comparison, 4D-Var outperforms the Kalman filters, achieving a DAB greater than 50% for observations with an error of 3 m.For longshore processes, due to its nature, the DAB is always greater than 80% because the shoreline position associated with longshore changes does not oscillate around an equilibrium position, as in the case of cross-shore processes, meaning that the background simulation completely diverges if no DA is applied.However, the PTWL is compromised in the cases of shorter pulse periods and noisy observations.In these cases, the best performing algorithm is the eKf, which achieves a PTWL greater than 80% for periods of 12 years or more with clean observations.For shorter periods or noisier observations, the parameter estimates from the eKf oscillate around the reference value but may fall outside the acceptable limits for the PTWL, as shown in figure 3(c).

Application to real cases
The different DA algorithms were tested at two coastal sites, Tairua and Castellón.Tairua is a crossshore-dominated 1.2 km long pocket beach located on the eastern coast of the Coromandel Peninsula, New Zealand.Shoreline, wave and water parameters were obtained from Montaño et al (2020).Shoreline observations were obtained from video cameras on a daily basis, with an acquisition error of 0.5 m.Castellón is a longshore-dominated beach located on the Mediterranean coast of Spain.Data were obtained from Álvarez-Cuesta et al (2021), and shoreline observations were extracted using the CoastSat algorithm on a bimonthly basis, with an acquisition error of 10 m (Vos et al 2019).
In this case, the no DA (background) and DA simulations were initialized using typical values for the model parameters following USACE (1984) and  Miller and Dean (2004) without assuming any prior system knowledge.The results are shown in figure 4. At Tairua, the DAB, calculated with respect to the background simulation during the validation period, exceeds 50% for 4D-Var and 44% and 36% for the eKf and SEnKf, respectively.All three algorithms provide satisfactory results, but 4D-Var outperforms the Kalman filters, yielding an RMSE of 4.4 m.At Castellón, after calibration with noisy satellite observations for longshore processes, 4D-Var is the worst performing algorithm, achieving a DAB of 71%, whereas the SEnKf reaches a DAB of 85%.The results are consistent with the conclusions derived from the synthetic cases for cross-shore and longshore processes.For cross-shore processes and almost clean observations, all the methods yield similar results, but 4D-Var slightly outperforms the Kalman filters.When modelling longshore processes using noisy observations, the Kalman filters perform better than 4D-Var, as the latter requires more accurate observations to achieve optimal performance.However, the DAB values obtained in the real cases are smaller than those obtained for the equivalent theoretical case shown in figure 2. This fact can be attributed to cascading forcing conditions and model errors arising from our incomplete knowledge of the physical system.

Discussion and conclusions
The results of the different DA methods are summarized in the form of a scorecard in table 1.For lowquality (noisy) data, the Kalman filters outperform 4D-Var for both cross-shore and longshore processes.
The strong constraint assumption of 4D-Var, which imposes the requirement that the system trajectory is determined by the model equations within the assimilation window, may need to be relaxed to properly integrate inaccurate data, as in weak-constraint 4D-Var (Kalnay et al 2007).However, 4D-Var and the eKf both provide better results than the SEnKf when the system is initialized far from the reference state and when it evolves nonlinearly for both cross-shore and longshore processes.In these cases, both 4D-Var and the eKf explicitly incorporate the tangent linear and adjoint operators of the system dynamics, while the SEnKf derives the corresponding magnitudes by propagating an ensemble of states.In such cases, the ensemble may not properly sample the true system covariance, resulting in filter divergence even if additive covariance inflation is applied to the ensemble members by considering the process noise covariance matrix Q.In contrast to the eKf and 4D-Var, the SEnKf is relatively easy to implement because this DA algorithm derives the required information from the physical system by propagating an ensemble of states without the need for further analytical derivations.This fact may explain why Kalman filters have thus far been favoured over variational methods in shoreline modelling.For simulations with reduced state vectors (O ∼10-100 elements), the eKf is more efficient than the SEnKf because it requires only one Ease of Implementation Cross-shore SEnKF forward simulation.For larger states, however, the propagation of the full analytically derived covariance matrix of the eKf is prohibitive, and its representation in the ensemble space of the SEnKf becomes a more affordable solution.
The results presented in this work seek to provide general guidelines to help decide which DA method to use on the basis of the main physical processes of interest as well as the available observations and system knowledge.While we have extensively examined conventional and widely used formulations for shoreline modelling, the adoption of alternative process-based or data-guided formulations would necessitate a comparable sensitivity analysis akin to the one presented here.This step is crucial for determining the most appropriate DA method for a given modelling approach.Additionally, similar conditions and simple rules have been used here to initialize all the methods in terms of the background error covariance B or process noise covariance Q.Although substantial work has been done regarding the initialization of B for state estimation problems (Bannister 2008a(Bannister , 2008b)), its definition for parameter estimation problems remains challenging (Smith et al 2013), as does the initialization of an ensemble consistent with prescribed B and Q matrices.The results may be substantially improved through careful fine tuning of those magnitudes, for instance, with a casespecific derivation of the matrix Q, as in Ibaceta et al (2020).Importantly, the use of DA does not guarantee a significant improvement in the results in real cases if the background covariance and process noise are not properly set.This is because the perfect model hypothesis is not satisfied and errors from wave dynamics and model parametrizations may cascade, hampering the convergence of the DA algorithms.Thus, further research is required to define the B and Q matrices and ensemble initialization for DA problems involving combined state and parameter estimation in shoreline modelling.
As remote sensing data are revolutionizing coastal science, this study aims to optimize the way such data are employed in combination with traditional physics-based equations to produce more accurate shoreline projections.Such optimization will lead to improved hazard forecasts for risk management and climate change adaptation.DA can overcome some of the main limitations of purely data-driven methods while also complementing them.The data dependence of data-driven methods can be compensated for by the consideration of the physics of the relevant processes in DA methods, while data-driven models can efficiently process raw observational data to feed into DA algorithms.DA is also an excellent tool for data handling, as it allows new observations to be directly integrated into real-time modelling systems, ultimately contributing to the development of digital twin models (Hoffmann et al 2023).

Figure 1 .
Figure 1.DA process following statistical and variational approaches.Please, note that the uncertainty in the state and in the observations are characterized by their covariance matrices.

Figure 2 .
Figure 2. Analysis of the role of observations on the performance of the DA algorithms for cross-shore and longshore processes.Panel (a): influence of the length of the observations campaign and the error of the observations.Panel (b): influence of the length of the observations campaign and the sampling interval.The DAB is shown in coloured squares and the inner dots represent the eK and CPU cost.

Figure 3 .
Figure 3. Analysis of the role of the system knowledge on the performance of the DA algorithms for cross-shore and longshore processes.Panel (a): influence of the prior knowledge of the system, measured as the length of the observations campaign and the initialisation quantile.Panel (b): influence of the error of the observations and the period of the non-stationary pulse of the true parameter's variation.Panel (c): time evolution of the shoreline position and the non-stationary evolution of the free parameters is displayed for cross-shore and longshore processes.The DAB is shown in coloured squares and the inner dots represent the eK in panel a and the PTWL in panel b.

Figure 4 .
Figure 4. Performance of the DA algorithms at two real sites: left, Tairua beach and right, Castellón beach.Top panels: the reference shoreline evolution is represented in black, the no DA simulation in grey, while in red, green and blue the DA shoreline predictions obtained using the SEnKf, the eKf and the 4D-Var, respectively.The bottom panels represent the temporal evolution of the parameters according to the different DA algorithms.The grey shaded area corresponds to the calibration period, while the reminder corresponds to the validation period.

Table 1 .
Score card of the different DA algorithms for cross-shore and longshore processes.The more stars, the better the performance of the algorithm.