Aurelien Callens
email: aurelien.callens@univ-pau.fr

Denis Morichon
email: denis.morichon@univ-pau.fr

Pedro Liria
email: pliria@azti.es

Irati Epelde
email: iepelde@azti.es

Benoit Liquet
email: benoit.liquet@univ-pau.fr

Magaly Koch

Yukiharu Hisaki

Xiaofeng Li

Zhixiang Fang

Quanyi Huang

Jaroslaw Tęgowski

Y-G Wang

L Fu

Automatic Creation of Storm Impact Database Based on Video Monitoring and Convolutional Neural Networks

Keywords: Coastal flooding, Coastal risk, Deep learning, Early warning system, Machine learning, Monitoring network, Statistical learning methods ANN Artificial neural networks. BN Bayesian networks. CNN Convolutional neural networks. DBN Dynamic bayesian networks. DRR Disaster risk reduction. EWS Early warning system. SLM Statistical learning methods. SNN Shallow neural networks. SOM Self organizing maps Artificial neural networks, Data assimilation, Error prediction, Gradient boosting trees, Random forest, Wave forecasting convolutional neural networks, storm impact database, transfer learning, video monitoring trained only on observational data Bayesian networks, Coastal flooding, Observational data, Statistical learning methods Heteroscedasticy, Model selection, Robust estimation, Temporal correlations

Over the last decades, the quantity of data related to coastal risk has greatly increased with the installation of numerous monitoring networks. In this era of big data, the use of statistical learning methods (SLM) in the development of local predictive models becomes more legitimate and justified. The objective of this thesis is to demonstrate how SLM can contribute to the improvement of coastal risk assessment tools and to the development of an early warning system which aims to reduce coastal flooding risk.

Three methodologies have been developed and tested on real study sites.

The first methodology aims to improve the local wave forecast made by spectral wave model with machine learning methods and data from monitoring networks. We showed that data assimilation with machine learning methods improve significantly the forecast of wave parameters especially the wave height and period. The second methodology concerns the creation of storm impact databases. Even though these databases are essential for the disaster risk reduction process they are rare and sparse. We therefore proposed a methodology based on a deep learning method (convolutional neural networks) to generate automatically qualitative storm impact data from images provided by video monitoring stations installed on the coast. The last methodology is about the development of a storm impact model with a statistical method (bayesian network) based exclusively on data acquired with diverse monitoring networks. With this methodology we were able to predict qualitatively the storm impact on our study site, the Grande Plage of Biarritz.

Résumé

Au cours des dernières décennies, la quantité de données relatives aux risques côtiers a fortement augmenté avec l'installation de nombreux réseaux de surveillance. Dans cette ère de big data, l'utilisation de méthodes d'apprentissage statistique dans le développement de modèles prédictifs locaux devient de plus en plus légitime et justifiée. L'objectif de cette thèse est de démontrer comment les méthodes d'apprentissage statistique peuvent contribuer à l'amélioration des outils d'évaluation des risques côtiers et au développement d'un système d'alerte précoce qui vise à réduire le risque d'inondation côtière.

Trois méthodologies ont été développées et testées sur différent sites d'étude. La première méthodologie vise à améliorer les prévisions locales de vagues faites par un modèle spectral de vagues avec des méthodes d'apprentissage automatique et des données provenant de réseaux de surveillance. Nous avons montré que l'assimilation de données avec des méthodes d'apprentissage automatique améliore de manière significative la prévision des paramètres des vagues, en particulier la hauteur et la période des vagues.

La deuxième méthodologie concerne la création de bases de données sur l'impact des tempêtes. Bien que ces bases de données soient essentielles dans le processus de réduction des risques de catastrophes, elles sont rares et peu nombreuses. Nous avons donc proposé une méthodologie basée sur une méthode d'apprentissage profond (réseaux de neurones convolutifs) pour générer automatiquement des données qualitatives sur l'impact des tempêtes à partir d'images fournies par des stations de surveillance vidéo installées sur les côtes. La dernière méthodologie concerne le développement d'un modèle d'impact de tempêtes avec une méthode statistique (réseau bayésien) basée exclusivement sur des données acquises avec divers réseaux de surveillance.

Grâce à cette méthodologie, nous avons pu prédire de manière qualitative l'impact de tempêtes sur notre site d'étude, la Grande Plage de Biarritz.

Mots-clés: Apprentissage automatique; Apprentissage profond; Inondations côtières; Méthodes d'apprentissage statistique; Réseau de surveillance; Risque côtier; Système d'alerte précoce.

First and foremost, I would like to express my sincere thanks and gratitude to my PhD supervisors, Benoit Liquet and Denis Morichon, for giving me their trust and proposing me numerous opportunities over the past years. I would like to thank Benoit in particular for his continuous cheerfulness and broad knowledge about statistics that made our weekly meetings always rich and more than pleasant.

I would also like to thank Denis who shared with me his deep knowledge about coastal engineering which was not my area of expertise and who helped me improve significantly my academic writing with his remarks and numerous proof-readings.

Their guidance is without a doubt a key ingredient of this thesis.

Beside my advisors, I would like to thank Deborah Idier and Pierre Lafaye de Micheaux for their constructive remarks and suggestions about my work and for taking the time to read and review this thesis.

I would like to express my gratitude to the persons who collaborated with me during this PhD thesis and who made it all the more interesting: Stéphane Abadie, Matthias Delpey, Bruno Castelle, Agnès Petrau, Maialen Sagarduy, Pedro Liria, Irati Epelde. I would like to give a special thank to all the personal of the doctoral school and the Laboratory of Mathematics and their applications of Pau, who helped me in any way for the administrative tasks and for the inscription of training courses.

I also thank my PhD and postdoc colleagues: Fania, Claire, Sebastien, Floren, Bastien, Téo for the stimulating discussions we had and the good time we shared during our meetings or during training courses.

From a more personal point of view, I am extremely grateful to my parents, my brother and my grandparents for their love and continuous support throughout my life and especially during this important period. I am also grateful to my friends Arthur, Charly for the gaming sessions spent together and their humor that always cheered me up. Finally, I am forever thankful to Mylène, my love and partner in life, for reassuring and supporting me with so much kindness and patience at every step of this project.

General context

Coastal areas have always been attractive habitats for human communities because of the abundant resources and valuable services provided by the diverse ecosystems associated with these areas [START_REF] Barbier | The value of estuarine and coastal ecosystem services[END_REF][START_REF] Mehvar | Quantifying economic value of coastal ecosystem services: a review[END_REF] and their strategic positioning at the interface between land and sea which offers access points to marine trade and transport. Over the last decades, the coastal zones have known an active urbanization [START_REF] Seto | A metaanalysis of global urban land expansion[END_REF] and population growth. Nowadays, the population density in coastal areas is significantly higher than in inland areas [START_REF] Neumann | Future coastal population growth and exposure to sea-level rise and coastal flooding-a global assessment[END_REF]. At world level, the coastal areas represent only 4% of the earth's total land area [START_REF] Barbier | Valuing ecosystem services for coastal wetland protection and restoration: Progress and challenges[END_REF], whereas they host 30% of the World's population [START_REF] Mea | Ecosystems and human well-being: wetlands and water synthesis[END_REF].

The same phenomenon is observed for Europe (Figure 1.1) where it is estimated that one third of the population live within 500 meters of the European seas or Oceans (European Commission, 2012).

The coastal zones also play a major role in the economy of human communities. As evoked earlier, their strategic positioning enables the construction of ports which are essential for international trade and transport. However, most of the economic value of these areas comes from the resources and services provided by the ecosystems present in these areas [START_REF] Barbier | The value of estuarine and coastal ecosystem services[END_REF]. The main valuable services are tourism, recreation and storm protection services [START_REF] Mehvar | Quantifying economic value of coastal ecosystem services: a review[END_REF]. In Europe, one third of the gross domestic product is produced within 500 meters of European seas and the economic value of these areas has been estimated between 500 and 1,000 billion Euros (European Commission, 2012).

Coastal risk

Coastal areas are exposed to different types of natural hazards such as storm waves, flood and erosion that can pose significant risks to assets and communities installed in these areas. According to the terminology of the United Nations relating to disaster risk reduction (UNGA, 2016), a risk is referred to the potential loss of life, injury, or destroyed or damaged assets which could occur to a system, society or a community in a specific period of time, determined probabilistically as a function of hazard, exposure, vulnerability (Figure In environmental studies, hazard usually refers to a phenomenon with given intensity and frequency that may cause loss of life, property damage or economic disruption (non exhaustive list). In the case of coastal areas, the hazards are mainly natural hazards such as storm waves, coastal flooding or erosion. Two categories of hazards can be distinguished: rapid onset hazard that occurs at a time scale of days to weeks and slow onset hazard that occurs at a time scale of decades to centuries [START_REF] Hill | Oceanography. Geosciences[END_REF]. Exposure alludes to the communities, the properties or other human assets which are subject to potential losses because of their location in hazard-prone areas [START_REF] Unga | Report of the open-ended intergovernmental expert working group on indicators and terminology relating to disaster risk reduction[END_REF]. Finally, vulnerability corresponds to the physical, social, economic and environmental factors associated with these exposed communities, properties or other human assets which increase their susceptibility to the impact of hazards.

Coastal flooding

In this thesis, a particular attention is paid to coastal flooding occurring on highly urbanized cities. This type of hazard falls in the class of rapid onset hazards. A flood is a general and temporary inundation of normally dry land areas. It is considered as coastal flooding when coastal processes such as tide, waves and storm surge are involved (Yang and Liu, 2020). This hazard occurs during extreme water-level events when the total water level exceeds the elevation of a defense infrastructure. It can cause severe damages on assets or buildings located behind the defense infrastructures (Figure 1.3). The total water level is a combination of several components, namely the astronomical tide, the storm surge and the wave runup (Figure 1.4).

Astronomical tide

The astronomical tide is one of the major components of the total water level. Tides are the rising and falling of the surface of the ocean caused by the combined gravitational attraction of the sun and the moon. In some regions, the tide level will determine whether or not there will be a flooding. For instance, it has been shown that coastal floodings occurred most of the times during high tides in the central part of the bay of Biscay [START_REF] Breilh | How frequent is storm-induced flooding in the central part of the Bay of Biscay?[END_REF].

Spring tides can also be an aggravating factor in coastal flooding. Twice a lunar month, the Earth, sun, and moon are nearly in alignment (full moon and new moon) which results in the so-called spring tides that are tides with slightly larger tidal range. A storm event occurring during a high spring tide is more likely to cause severe floods.

Storm surge

A sea surge corresponds to the difference between the observed sea level and the astronomical tidal level [START_REF] Pirazzoli | Surges, atmospheric pressure and wind change and flooding probability on the atlantic coast of france[END_REF]. This difference is the result of the interactions between the sea and several forcings such as the atmospheric pressure, the winds or the waves. Atmospheric pressure has a significant effect on the surge level. This phenomenon is known as inverse barometric effect.

It has been estimated that a decrease of 1 hPa in atmospheric pressure leads to an elevation of 1 cm in total water level (Harris, 1963). Winds can also move water masses toward the coast and lead to increased surge level. In shallow coastal areas, winds have generally stronger effect on sea surge than atmospheric pressure (Arnaud and Bertin, 2014). Finally, it has been proven that sea state can influence significantly the surge level (Bertin et al., 2015).

During a storm, the extreme conditions (low pressure, high winds and large waves) lead to a significant surge which is qualified as storm surge.

Wind waves induced components

By breaking on the shore, waves have a direct effect on the total water level. This effect is named "wave runup" and is defined as the maximum vertical extent of wave uprush on a beach or structure above the still water level. The wave runup is the sum of two components: "wave setup" which is a nearly static component and "swash" which is a dynamic component at the wave scale. Wave setup is the increase in mean water level above the stillwater level due to momentum transfer to the water column by waves that are breaking or otherwise dissipating their energy. The swash corresponds to the vertical elevation of the lens of water that washes up on the beach after a wave has broken. The swash dynamics is controlled by incident waves (period smaller than 20s) and infragravity waves (period ranging between 30 and 300s) (Stockdon et al., 2006). The contribution of each component varies with the beach type in such a way that swash in dissipative beaches is dominated by infragravity waves, and in reflective beaches by incident waves [START_REF] Plant | How well can wave runup be predicted? comment on laudier et al[END_REF].

During a storm event, wave height and period are generally larger than average. This can cause a significant elevation of the total water level because wave height and period are positively correlated with wave runup [START_REF] Diwedar | Investigating the effect of wave parameters on wave runup[END_REF]. It has been shown that the contribution of waves to the total water level is more important than the contribution of winds in zones where the continental shelf is reduced [START_REF] Kennedy | Tropical cyclone inundation potential on the Hawaiian Islands of Oahu and Kauai[END_REF]. Waves can also interact with other components of the total water level. Several authors have demonstrated that waves can influence the storm surge [START_REF] Nicolle | Amplification of the storm surges in shallow waters of the pertuis charentais (bay of biscay, france)[END_REF][START_REF] Bertin | Importance of wave age and resonance in storm surges: The case Xynthia, Bay of Biscay[END_REF]Bertin et al., , 2015)). Indeed, the sea state can enhance the wind friction, inducing the movement of larger bodies of water toward the coast which results in larger storm surge.

Coastal flooding and global change

The risk related to coastal flooding is expected to increase significantly in the future. On one hand, global warming will have an impact on all the components of the total water level which are responsible for coastal flooding.

Since the beginning of the last century, the global mean sea level has been rising at an accelerated rate [START_REF] Chen | The increasing rate of global mean sea-level rise during 1993-2014[END_REF] and is expected to keep rising during the next decades, reaching an increase of 50 cm (RCP 4.5) to 80 cm (RCP 8.5) by the end of the century [START_REF] Kopp | Probabilistic 21st and 22nd century sea-level projections at a global network of tide-gauge sites[END_REF][START_REF] Mengel | Future sea level rise constrained by observations and long-term commitment[END_REF].

The sea level rise and the global warming will induce changes in the amplitudes and phases of tides [START_REF] Idier | Sea-level rise impacts on the tides of the European Shelf[END_REF][START_REF] Pickering | The impact of future sea-level rise on the global tides[END_REF]. Even though an increase in storminess in the north Atlantic for the future decades is still debatable [START_REF] Bengtsson | Will extratropical storms intensify in a warmer climate[END_REF][START_REF] Feser | Storminess over the North Atlantic and northwestern Europe-A review[END_REF], wind waves [START_REF] Perez | Statistical multi-model climate projections of surface ocean waves in Europe[END_REF] and storm surges [START_REF] Marcos | Changes in storm surges in southern europe from a regional model under climate change scenarios[END_REF][START_REF] Little | Joint projections of us east coast sea level and storm surge[END_REF] will be affected by global warming. These changes, especially in mean sea level, will result in more frequent and intense coastal floodings [START_REF] Vitousek | Doubling of coastal flooding frequency within decades due to sea-level rise[END_REF][START_REF] Vousdoukas | Global probabilistic projections of extreme sea levels show intensification of coastal flood hazard[END_REF]Taherkhani et al., 2020).

On the other hand, the stakes on coastal areas will keep growing. More people will be exposed to coastal flooding as the populations in these zones are expected to increase significantly in the future decades [START_REF] Neumann | Future coastal population growth and exposure to sea-level rise and coastal flooding-a global assessment[END_REF]. Urban areas will keep developing in coastal zones to welcome the growing populations resulting in more assets, buildings and infrastructures to protect from coastal flooding. In the study of Vousdoukas et al. (2018a), it has been esti-mated that the current expected annual number of people exposed to coastal flooding of 102 000 will reach 1.52-3.65 million by the end of the century for Europe depending on the scenario considered. Concerning the current expected annual damage for Europe, the authors found an increase by two to three orders of magnitude by the end of the century if nothing is done (from 1.25 billion to 93-961 billion euros depending on the scenario used).

Reducing coastal risk

Given the current flooding risk and its future increase in coastal areas, efficient disaster management is needed to identify, assess or reduce the risks of disaster.

Generally, the process of Disaster risk reduction (DRR) is represented as a cycle (Figure 1.5) with 4 steps: response, recovery, prevention, preparedness. Yang and Liu (2020) The response corresponds to the actions and measures taken directly before, during or immediately after a disaster in order to save lives, reduce health impacts, ensure public safety (UNGA, 2016). The recovery is the step occurring in the post-disaster phase, and corresponds to the restoring of the economic, social, environmental assets (non exhaustive list) affected by this hazard. Prevention and preparedness steps occur during the pre-disaster phase.

The prevention step corresponds to all the activities and measures taken to avoid existing and new disaster risks. Finally, the preparedness step represents the knowledge and capacities developed by governments, communities to effectively anticipate, respond to the impacts of likely, imminent or current disasters [START_REF] Unga | Report of the open-ended intergovernmental expert working group on indicators and terminology relating to disaster risk reduction[END_REF]. Among all the measures employed during the disaster management process, we can distinguish two families with different aims: mitigation or preparedness.

The measures aiming to mitigate directly the impact of a coastal hazard are often structural measures that correspond to engineering constructions. The measures aiming to prepare and prevent a disaster are often nonstructural measures such as actions, legislation or warning systems [START_REF] Meyer | Economic evaluation of structural and non-structural flood risk management measures: examples from the mulde river[END_REF].

Our interest lies in these measures as they are more proactive, relatively cost effective and require a short time to implement.

Preparedness

The measures aiming for preparedness focus on two key factors which are the knowledge about the hazards and the ability to predict in advance such events to give an early warning to the communities.

Coastal monitoring networks

The most efficient way to build knowledge about hazards is to install permanent monitoring systems. Over the last decades, numerous coastal monitoring systems have been installed along the coasts in order to study the processes that are responsible for coastal flooding. These systems are generally regrouped into networks to provide continuous and sustainable data.

Along the french coasts, we find the RONIN network operated by SHOM that regroups tidal gauges which are measuring the total water level, the CAN-DHIS network operated by CEREMA constituted by numerous offshore directional wave buoys which record the wave characteristics and a network of weather stations operated by Metéo-France. There are also remote sensing systems such as video monitoring station or satellites. Video monitoring systems are installed on the beach to study coastal processes such as beach morphology changes, wave runup, and coastal currents (Splinter et al., 2018;Buscombe and Carini, 2019). Examples of networks for video monitoring systems are euskoos (https://www.kostasystem.com/fr/) for the Basque coast, WebCAT (https://secoora.org/web-cameras/) for the east coast of USA or the coastal imaging network of the water research laboratory of New South Wales (http://ci.wrl.unsw.edu.au/) for australian coasts. Satellites are employed to perform various tasks such as mapping ocean currents [START_REF] Klemas | Remote sensing of coastal and ocean currents: An overview[END_REF], estimating wave parameters [START_REF] Shao | Ocean wave parameters retrieval from sentinel-1 sar imagery[END_REF][START_REF] He | Remote sensing of ocean waves by polarimetric sar[END_REF] and also monitoring meteorological conditions [START_REF] Kidd | A review of satellite meteorology and climatology at the start of the twenty-first century[END_REF]. During disaster or post-disaster phase, remote systems can assess the impact of the hazard which can be precious for building Early warning system (EWS) [START_REF] Klemas | The role of remote sensing in predicting and determining coastal storm impacts[END_REF].

Data collected from monitoring systems are essential for the preparedness as they are used as input for statistical learning models or EWS (Valchev et al., 2014;Van Dongeren et al., 2018). They are also used as ground truth to calibrate and validate numerical wind wave models [START_REF] Lavidas | Application of numerical wave models at european coastlines: A review[END_REF].

Predictive models and Early warning systems

To complement coastal monitoring systems, the development of accurate predictive models and Early warning system (EWS) are crucial measures for the preparedness step. According to [START_REF] Unga | Report of the open-ended intergovernmental expert working group on indicators and terminology relating to disaster risk reduction[END_REF], an EWS designate a system including four interrelated key elements: (1) disaster risk knowledge based on the systematic collection of data and disaster risk assessments; (2) monitoring, analysis and forecasting of the hazards and possible consequences; (3) dissemination and communication of timely accurate warnings and associated information on likelihood and impact; and ( 4) preparedness at all levels to respond to the warnings received. The development of EWS is considered as the most cost-effective measure, as EWS save both lives and properties [START_REF] Rogers | Costs and benefits of early warning systems[END_REF].

Various methods can be used inside an EWS to forecast the characteris-tics of a hazard or predict its impact. These methods can be divided into 2 main groups: physics-based methods (numerical modeling) and data-driven methods (statistical learning). Physic-based models are representations of physical processes in mathematical terms. They are usually composed by one or more governing equations that are based on theory and fundamental knowledge. Because these equations are often impossible to solve or because initial or boundary conditions are not known [START_REF] Larson | Numerical Modeling[END_REF], numerical modeling techniques are usually employed to approximate solutions. On the opposite, data-driven approaches are based on the analysis of the data about a specific system. They aim to find relationships between the system variables (input and output) without explicit knowledge of the physical behavior of the system [START_REF] Solomatine | Data-driven modelling: some past experiences and new approaches[END_REF].

This thesis focuses on the data-driven methods which are getting more attention over the last decade because of the increasing amount of data collected by the monitoring networks and the rapid development of data science.

Statistical learning methods

Statistical learning is one of the data-driven approaches. It is a field combining machine learning and statistics, it regroups tools and methods for modeling, predicting and understanding complex data. SLM aim to solve different type of learning problems such as unsupervised learning which aims to find the structure of the data or supervised learning which aims to map data to a desired output. The latter is our main focus because it has applications in predictive modeling. Over the past year, the field of statistical learning has gained increasing interest and attention as can be seen from the Calculation of the normalized measure is detailed here: https://support. google.com/trends/answer/4365533?hl=en.

Our interest will lie upon the use of three categories of SLM namely machine learning, deep learning and statistical methods in coastal risk analysis.

Machine learning methods are algorithms that learn pattern from pre-defined data features in order to make the most accurate predictions for new data.

The features extraction in machine learning is a manual process that requires domain knowledge. Although they are considered as black box models, they are widely employed due to their prediction accuracy [START_REF] Burkart | A survey on the explainability of supervised machine learning[END_REF]. Deep learning is sub-class of machine learning methods based on deep artificial neural networks which are algorithms inspired by the structure and function of the brain. Deep learning methods aim to discover the mapping from data features to the target variable but also discover which features to extract from data contrary to machine learning methods. Deep learning methods generally require a large amount of data and computational power to train.

Finally, statistical methods are methods based on statistic, they can be used as predictive models but they are more oriented toward the understanding of the data and the study of relationships between variables.

Machine learning

In coastal risk analysis, machine learning methods are commonly employed due to their predictive ability. They are employed to perform short-term forecast of coastal processes responsible for coastal flooding or to improve forecast of the same processes made by process-based numerical models.

Concerning these two applications, we observed in the literature a preference for Shallow neural networks (SNN) as a machine learning algorithm.

Indeed, they are used as off-the-self method in numerous works. SNN designate a class of feed-forward Artificial neural networks (ANN). They are generally simple ANN with three layers of nodes: an input layer, one hidden layer and an output layer. They have been employed to obtain short-term forecast of wave parameters, storm surge and tide level (Table 1.1). The explanatory variables (input) are mainly measurements obtained from monitoring networks (buoy, weather stations, tide gauge) or predictions made by numerical wave models.

They have also been employed to improve the predictions made by numerical wind wave models with a data assimilation procedure called error prediction method (Babovic et al., 2001;Makarynskyy et al., 2005;Moeini et al., 2012;Deshmukh et al., 2016;Londhe et al., 2016). In this case, the explanatory variables are wave parameters simulated by the wind wave model and meteorological variables and the response variable are the deviations of the wind wave model. In Londhe et al. (2016), this data assimilation method with SNN led to a significant improvement of significant wave height forecast, with the correction of the underestimation (no remaining bias after the assimilation) and RMSE in average 50 % lower.

Table 1.1: Scientific works using SNN in the predictive modeling of tide, waves and storm surge.

Modeled parameter References

Tide [START_REF] Lee | Back-propagation neural network for long-term tidal predictions[END_REF], [START_REF] Makarynskyy | Predicting sea level variations with artificial neural networks at Hillarys Boat Harbour, Western Australia[END_REF], [START_REF] Makarynska | Predicting sea-level variations at the Cocos (Keeling) Islands with artificial neural networks[END_REF], [START_REF] Granata | Artificial intelligence models for prediction of the tide level in venice[END_REF] Wave height Deo et al. (2001), [START_REF] Tsai | Neural network for wave forecasting among multi-stations[END_REF], Makarynskyy (2005), Mandal and Prabaharan (2006), [START_REF] Browne | Near-shore swell estimation from a global wind-wave model: Spectral process, linear, and artificial neural network models[END_REF] Storm surge level [START_REF] Lee | Neural network prediction of a storm surge[END_REF], [START_REF] Tseng | Application of artificial neural networks in typhoon surge forecasting[END_REF], [START_REF] You | Storm surge prediction using an artificial neural network model and cluster analysis[END_REF], [START_REF] Hashemi | An efficient artificial intelligence model for prediction of tropical storm surge[END_REF], [START_REF] Bezuglov | Multi-Output Artificial Neural Network for Storm Surge Prediction in North Carolina[END_REF], [START_REF] Kim | A Surrogate Modeling for Storm Surge Prediction Using an Artificial Neural Network[END_REF], [START_REF] Lee | The Study for Storm Surge Prediction Using Generalized Regression Neural Networks[END_REF], [START_REF] Quintana | Statistical forecast of the marine surge[END_REF] Despite the popularity of SNN, other machine learning methods have been employed to predict coastal processes responsible for coastal flooding. In the literature, significant wave height has been predicted by support vector machines [START_REF] Elgohary | Significant deep wave height prediction by using support vector machine approach (alexandria as case of study)[END_REF][START_REF] Berbić | Application of neural networks and support vector machine for significant wave height prediction[END_REF], regression trees (Etemad-Shahidi and Mahjoobi, 2009) and random forest [START_REF] Mafi | Forecasting hurricane wave height in gulf of mexico using soft computing methods[END_REF].

The storm surge has been also predicted by different algorithms such as support vector machines [START_REF] Hashemi | An efficient artificial intelligence model for prediction of tropical storm surge[END_REF], regression trees and random forest [START_REF] Granata | Artificial intelligence models for prediction of the tide level in venice[END_REF]. Over the last decades, an increasing number of works have been comparing the SNN with other machine learning methods proving that this method is not the best for all problems. In [START_REF] Granata | Artificial intelligence models for prediction of the tide level in venice[END_REF], they showed that regression trees outperforms SNN in most of the case for the tide forecast. Same observation was made on the prediction of significant wave height [START_REF] Etemad-Shahidi | Comparison between m5 model tree and neural networks for prediction of significant wave height in lake superior[END_REF].

Deep learning

With the increasing amount of data and computational power, deep learning methods such as recurrent neural networks have been increasingly used over the last years in the predictive modeling of coastal processes. Recurrent neural networks are deep neural networks specialized in the sequences and time series prediction (Bengio et al., 2017). This type of networks have been used successfully to accurately predict time series of significant wave height (Mandal and Prabaharan, 2006;[START_REF] Sadeghifar | Coastal wave height prediction using recurrent neural networks (rnns) in the south caspian sea[END_REF][START_REF] Savitha | Regional ocean wave height prediction using sequential learning neural networks[END_REF][START_REF] Alqushaibi | Enhanced weight-optimized recurrent neural networks based on sine cosine algorithm for wave height prediction[END_REF] or to predict storm surge level [START_REF] Di Nunno | Forecasting of extreme storm tide events using narx neural network-based models[END_REF][START_REF] Quintana | Statistical forecast of the marine surge[END_REF]. More recently, they have been used to improve the predictions of a numerical wind wave model by integrating both the local data and the temporality in the errors of the numerical wind wave model [START_REF] Zhang | Ocean wave height series prediction with numerical long short-term memory[END_REF].

Convolutional neural networks (CNN) which are deep neural networks specialized in images analysis are also more and more employed to analyze the large quantity of images created by the video monitoring stations located on the shore. However, there are only few applications of CNN based on the images of the video monitoring stations such as the estimation of the nearshore bathymetry [START_REF] Benshila | A deep learning approach for estimation of the nearshore bathymetry[END_REF] or wave-tracking [START_REF] Kim | Wave-tracking in the surf zone using coastal video imagery with deep neural networks[END_REF]. So far, they do not have a direct application for the analysis of coastal risk.

Statistical methods

In the study of coastal risks, statistical methods are used for their ability to bring understanding about the data and to study the relationships between variables. Among the statistical methods, Bayesian networks (BN), a class of probabilistic graphical models, are appreciated methods in coastal risk modeling. Their low computational cost, their ability to represent complex systems by integrating different sources of data and their intuitive representation are non negligible advantages compared to other modeling approaches such as numerical process based models.

Most of the applications of BN in coastal engineering domain concern the translation of forcing variables (e.g. wave, weather and tide conditions) into impact and damages on the shore during storm events. For instance, they have been employed to predict coastal cliff erosion (Hapke and Plant, 2010), shoreline retreat (Beuzen et al., 2018), dune retreat and erosion (Palmsten et al., 2014;den Heijer et al., 2012) and barrier island response (Plant and Stockdon, 2012;Wilson et al., 2015) resulting from coastal storms. They are also used as surrogates of process-based models. In Poelhekke et al. (2016) and Plomaritis et al. (2018), BN are trained using output data from many precomputed process-based model simulations. Once trained, the networks can be conditioned with forecast of the hydraulic boundary conditions to obtain instantaneously forecast of onshore hazards. By avoiding the computation time associated with process based models, bayesian networks are a great assets for operational EWS.

Advantages and limitations of SLM

The main advantage of SLM over process-based models is the computation time. The training time for such methods is usually in the order of minutes (or hours for deep learning) and the predictions made by these methods are instantaneous compared to numerical process-based models. In addition, SLM are relatively simple to use and with enough data they can achieve similar or even higher performances than numerical wind wave models for site-specific wave parameters [START_REF] Browne | Near-shore swell estimation from a global wind-wave model: Spectral process, linear, and artificial neural network models[END_REF][START_REF] Savitha | Regional ocean wave height prediction using sequential learning neural networks[END_REF].

The main limitation of these methods is that their performances depend directly on the quantity and quality of the data. Deployment of these models are only possible in sites equipped with monitoring systems. In the future, data availability will not be as much as a problem. Indeed, monitoring networks are expected to be more dense as new technological advances will improve sensors and their connectivity, and as the cost of the components will decrease [START_REF] Marcelli | Toward the widespread application of low-cost technologies in coastal ocean observing (internet of things for the ocean)[END_REF]. In this future context, the use of statistical methods adapted to big data will be necessary [START_REF] Xu | Internet of things in marine environment monitoring: A review[END_REF].

Problematic and objectives of the thesis

Considering that coastal flooding risk is expected to increase in the next decades and that SLM offers a promising potential in an era of big data, this thesis aims to answer the following question:

How SLM can contribute to the improvement of coastal risk assessment tools and to the development of an EWS which aims to reduce coastal flooding risk ?

This problematic is broken down into 3 sub-questions that will be investigated in this thesis:

• How to improve, with SLM, the local forecast of spectral wave models that are known to underestimate wave parameters during extreme events?

• Can we employ SLM to constitute storm impact databases that are rare and not routinely collected but are essential in the disaster risk reduction process?

• Can we fully develop an accurate predictive model linking offshore hydraulic boundary conditions into onshore hazards based on data collected from monitoring network and SLM?

Outline of the thesis

The presentation of the research carried out in this thesis is organized in 4

chapters each dedicated to a specific scientific issue linked to the application of SLM to the study of coastal flooding and its prediction (Figure 1. 7).

Chapter 1

Gain local understanding with SLM

Chapter 2

Improve wave forecast with SLM and local observations

Chapter 3

Creation of storm impact database based on video monitoring and deep learning

Chapter 4

Bayesian networks to model storm Impact based on observational data and SLM

Local knowledge and Supervised models

Data assimilation on wave data Storm impact database Figure 1.7: Connection between the chapters of this thesis.

The first chapter acts as an introduction on SLM and their applications in the study of coastal processes. The objective of this chapter is to demonstrate the interest of using these methods by highlighting the fact that they can provide knowledge about a coastal process at a local scale. Local knowledge is essential in the development of an EWS. This statement is especially true if the EWS is built with a statistical learning method such as bayesian networks, which is the case in this thesis. This chapter is organized in three parts, each presenting a different application of SLM in the study of coastal processes. Two of them concern supervised learning with the modeling of two coastal processes related to coastal flooding: storm surge level and wave runup on the beach. The last application of SLM concerns unsupervised learning and shows how to characterize local wave climate with clustering method.

In the second chapter of this thesis, we demonstrate the ability of SLM to improve local wave predictions made by numerical wind wave models. The accuracy of the data used in the development of the EWS is crucial to have reliable warnings. However, it has been shown that numerical wind wave models have a tendency to underestimate certain wave parameters during stormy conditions. We present alternatives SLM (random forest and gradient boosting trees) to neural networks which are traditionally used in the literature to perform data assimilation. We demonstrate that these methods, that have never been used for this task, leads to better improvements than neural networks in addition to providing explicability with the predictive power of the variables.

This chapter is constituted by an article published in Applied Ocean Research.

The third chapter demonstrates that SLM, especially deep learning methods, can be used to create automatically a storm impact database with images from video monitoring networks. Building a bayesian network that can be employed in EWS to predict coastal flooding requires a lot of data. It needs data about the coastal processes that are responsible for coastal flooding (tide, waves, storm surge, meteorological conditions) but it also needs data about the impact of coastal flooding. Unlike data about the coastal processes, data about the impact of coastal flooding are rare and sparse. Convolutional neural networks, which are deep learning methods, are used to classify the video monitoring images into three storm impact regimes which are categories of coastal flooding risk. Once trained, these networks can predict the storm impact regimes of newly created timestacks, generating an incremental storm impact database. This chapter is composed by an article published in Remote Sensing.

The fourth chapter shows that SLM can be employed to predict coastal flooding risk. This chapter regroups the models and data acquired through the previous chapters and assemble them in a Bayesian network that can be employed in a EWS to predict coastal flooding risk (Figure 1. 7). In this chapter, we evaluate the predictive performance of a BN exclusively based on observational data from diverse monitoring networks. In addition, we propose a methodology based on SLM to extend the storm impact and atmospheric surge data that are limited for the study site of Biarritz. A second BN is trained on the extended database and its performances are compared with the first BN to see the gain of performances associated with the extension of the database. This chapter is constituted by an article in preparation for a submission to Natural Hazards.

Finally, we conclude this thesis by answering the main research question, summarizing the scientific contributions associated with this thesis and discussing the potential perspectives. The second part of this chapter is about the modeling of wave runup on the beach depending on the characteristics of the offshore waves. This part is quite similar to the first part of chapter one in terms of methodology, however, an additional statistical method is tested which is the generalized additive model. This method has the advantage to better illustrate the relationships between the variable of interest and the covariates (variables used to make the predictions).

The last part of this chapter aims to characterize the wave climate at Anglet's buoy location by studying the relationships between the local wave characteristics observed at the buoy and the synoptic weather (sea level pressure) over the North Atlantic ocean. Contrary to the two other parts of this chapter, the characterization of the wave climate is an unsupervised task and other statistical learning methods must be employed such as dimensionality reduction (PCA) or clustering methods (Self Organizing Maps).

Sea surge modeling

Storm surges play a major role during coastal floodings, it is therefore essential for coastal researchers and stakeholders to have predictive models for this phenomenon. Numerical modeling is commonly employed for storm surge modeling. For a general review on this subject the reader is referred to the review of [START_REF] Flather | Existing operational oceanography[END_REF] and for examples of operational models to the works of [START_REF] Lionello | A data assimilation procedure for operational prediction of storm surge in the northern Adriatic Sea[END_REF] or [START_REF] Souza | Application of an operational storm coastal impact forecasting system[END_REF]. Over the last few years, statistical learning models (especially shallow neural networks) have been employed to model the storm surge [START_REF] Bezuglov | Multi-Output Artificial Neural Network for Storm Surge Prediction in North Carolina[END_REF][START_REF] Hashemi | An efficient artificial intelligence model for prediction of tropical storm surge[END_REF][START_REF] Lee | The Study for Storm Surge Prediction Using Generalized Regression Neural Networks[END_REF]. These models have the advantage to be less computationally expensive than process-based numerical models but they requires large observational dataset to be trained.

In this first part, we aim to train a statistical learning model to predict precisely sea surge level at Socoa tide gauge (South west of France), while learning about the local phenomenon by studying variable importance. Variable importance is a measure of the predictive power of a variable in a model.

It can be used to sort variables from most to least predictive, allowing one to have more insights on the problem and to perform feature selection when there are too many variables. This model can be potentially used to predict operationally local storm surge and the insights brought by the variable importance analysis can be included in the development of an EWS based on Bayesian networks (BN). Most authors take into account the temporal effect for the modeling of the atmospheric surge level and therefore use SNN as the "off-the-shelf" method.

Data

In this work, we do not take into account the temporal effects. It is why we are comparing several methods to find the most appropriate one for our problem.

We tested usual supervised learning algorithms:

• Linear model

• Shallow neural networks • Random forests • Gradient boosting trees
The reader is referred to [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF] for details on neural networks, random forest and gradient boosting trees. The algorithms all have the same explanatory variables to predict the storm surge level: the astronomical tide level (harmonic analysis), the local meteorological conditions (average wind speed and direction above 10 meters, atmospheric pressure) and the sea state characteristics (significant wave height, peak wave period and wave direction).

The best statistical learning model is the one which obtains the best performances (lowest root mean squared errors) on the test set.

Bayesian optimization for hyperparameter tuning

For complex learning algorithms such as random forest, gradient boosting trees and neural networks, hyperparameter values must be optimized. Hyperparameters are parameters whose values are specified by the user before the training process begins, they affect the structure of a learning algorithm and how well it trains. They have a non negligible impact on the final results.

In this work, Bayesian optimization is employed to select the optimal hyperparameter values. Bayesian optimization is an iterative algorithm that aims to minimize an objective function, in our case the root mean squared errors (RMSE). First, it builds a probability model (Gaussian process) of the objective function. Then it uses this surrogate model to select the most promising values of hyperparameters to evaluate. Once the promising combination of values have been evaluated, the probability model is updated and searched again for the most promising combination. This process is repeated several times. This method is employed because it is very efficient for tuning hyperparameter values and it usually requires less iterations than other methods such as grid or random search (Bergstra et al., 2011). In depth details of this method are given in the works of Snoek et al. (2012); Marchant and Ramos (2012) and Shahriari et al. (2015).

To tune the hyperparameter values of our models, bayesian optimization is coupled with a 5-fold cross validation on the the training set. The objective function to minimize for the Bayesian optimization method is the average outof-sample RMSE. The Bayesian optimization for our data is performed using the R package Tidymodels. This package allows for the training and tuning of several models simultaneously. To find optimal hyperparameters values for the different models (random forest, gradient boosting trees and neural networks), random combinations of hyperparameter values are first evaluated to serve as search base for the bayesian method (5 in this study), then an acquisition function (upper confidence bound) is used to find the next combination values to evaluate (this step is repeated 15 times).

Permutation Feature Importance

There is a lot of methods to compute variable importance depending on the type of model used. A complete review about these methods is presented in the work of [START_REF] Wei | Variable importance analysis: A comprehensive review[END_REF]. Permutation Feature Importance is one of the methods used to measure variable importance and is the one used in this work.

This method was first described for random forest by [START_REF] Breiman | Random forests[END_REF] and was later adapted to other models by [START_REF] Fisher | All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously[END_REF].

This method consists in shuffling the values of a variable (process called permutation) and measuring the increase in the model's prediction errors associated with these permutations. A variable is considered as important if the permutations cause an increase in the model errors. This increase of model errors means that the model relied strongly on this given variable. On the contrary, a variable is considered as an unimportant variable if the permutations do not increase the model errors.

Permutation feature importance is implemented in numerous R package.

This method is usually implemented by default in the packages specialized for random forest and gradient boosting trees. The model agnostic version of [START_REF] Fisher | All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously[END_REF] used in this work comes from the R package vip.

Results

The The importance of the variables for random forest is displayed in figure 2.3.

The variable with highest predictive power is the atmospheric pressure. This is expected as the atmospheric pressure influences significantly the surge level (Harris, 1963). The second, third and fourth most important variables are related to wave characteristics. In coastal areas with a reduced continental shelf (such as our site), it has been proven that the wind effect is limited (hence its low predictive power) and the contribution of both the atmospheric pressure and the waves is higher [START_REF] Kennedy | Tropical cyclone inundation potential on the Hawaiian Islands of Oahu and Kauai[END_REF][START_REF] Chaumillon | Storm-induced marine flooding: Lessons from a multidisciplinary approach[END_REF].

Discussion

The main limitation of the SLM is the quantity of data. Training models with more data results in models that are able to generalize better from a higher amount of information and therefore results in better performances on unseen data. This is especially true as we are modeling sea surge which can become With its promising performances, this model could be used in operational settings to forecast sea surge by using wave and meteorological forecast as explanatory variables. However, this is not recommended as this model has been trained with observational data and not on forecast data that might contain potential bias. The best solution for an operational forecast model should be to train the model with wave and meteorological forecast data in order to take into account the potential bias of the forecast models.

In the future, this work could be extended by adding explanatory variables containing temporal aspect. This could be the values of atmospheric pressure at previous time steps. In this framework, artificial neural networks could perform better as they are known to handle efficiently time series. Other input variables could be also used to improve the modeling such as the meteorological or wave data at different locations near the site.

Conclusion

We trained statistical learning models to predict the sea surge for the site of Socoa using several explanatory variables including meteorological conditions, sea state characteristics and tide level. By comparing the performances of these models on the test set, we defined random forest as the best model for this site with RMSE value of 5.73 cm.

The importance of the variables for random forest were computed with the permutation method. The variable importance analysis gave us insights about the local sea surge and the order of importance corresponded to the characteristics of the site and to what was found in the literature in similar sites. The knowledge acquired during this analysis can be helpful in the development of an EWS based on BN.

Runup elevation on the beach

Coastal floodings originate from the interaction of several processes which all influence the total water level. Among these processes, there are the tide, the meteorological conditions and also the waves. The contribution of the waves to the total water level is called wave runup and corresponds to the maximum onshore elevation reached by a wave, relative to the wave-averaged shoreline position (Flanders Marine Institute, 2021). Wave characteristics are available through observations (buoy) or numerical modeling. However, the wave characteristics only gives a broad idea of the runup on the beach. To transform offshore wave characteristics into runup elevation on the beach, two modeling approaches can be employed: numerical modeling or data-driven approaches.

There are various numerical models that solve hydrodynamic equations to estimate runup elevation on the beach. More details on these different models can be found in the work of [START_REF] Fiedler | Numerical modeling of wave runup on steep and mildly sloping natural beaches[END_REF]. Runup on the beach can also be modeled using empirical formula or statistical models estimated with local data. Many empirical formulas have been developed over the year, a complete review is given in the work of da [START_REF] Da Silva | On the prediction of runup, setup and swash on beaches[END_REF]. A non negligible advantage for the data-driven methods is in the computational aspect, where the computational effort and time are lower than those of numerical modeling methods.

Measurements of the wave runup were performed on the study site of Biarritz (South west of France). By comparing the estimations obtained by the formula of Stockdon et al. (2006) (the most used empirical formula) with the measurements, large estimation errors were found, highlighting the fact that this empirical formula was not adapted to this site. An empirical formula based on a linear model was developed by the SIAME researchers. This second formula yielded runup estimation much closer to the observations (Figure 2.4). In this second part of the first chapter, we aim to extend the work made the Grande Plage of Biarritz. We want to train a statistical learning model that yields better runup predictions than the empirical formula developed on the Grande Plage de Biarritz, while learning about the local phenomenon.

Study site

The Grande Plage de Biarritz (GPB) is an urbanized beach located in the South Ouest of France. This beach has a high socio-economic value due to its location near the city center, its historical heritage and its tourist appeal (Morichon et al., 2018)

Data

Runup measurements

In the framework of the MAREA project, runup was measured on the Grande Plage de Biarritz by the SIAME (Applied sciences in mechanics and electrical engineering) laboratory. The measurement of wave runup were obtained using the video monitoring station installed on the Grande Plage de Biarritz. This video monitoring station consists in 4 cameras, controlled by the open source software SIRENA (Nieto et al., 2010). The pixel intensities were sampled along two transects (Figure 2.5) at 1Hz during 14 min to produce timestack images. Standard transformations from image to world coordinates (image orthorectification) were employed. Then, Otsu's thresholding method was applied to the collected timestack images to obtain time series of the cross-shore position of the waterline. This semi-automatic extraction method has already been employed to extract waterline in timestack images (Vousdoukas et al., 2012) and aims to separate pixels into two classes based on their intensities. The obtained cross-shore coordinates were transformed into elevations using the topographic information collected with a NRTK GNSS at low tide during the 14 days of field surveys. These 14 days of fields survey are spread between the 4th December 2017 and the 24th January 2019. During these 14 days, runup R and the slope of the beach (β) were measured every 15 minutes.

Wave parameters

The joining the runup measurements and the waves characteristics we obtain a dataset of 220 observations. Due to the small number of observations and the temporal discontinuity of the data, we do not consider the temporal aspect in the modeling step presented in the next section.

Methodology

Statistical learning methods

We compare several statistical learning methods with different characteristics: linear models, generalized additive models (GAM) and random forest.

Linear models are the simplest method and will serve to determine baseline performance. Random forest are ensemble learning methods which rely on weak learners (classification and regression trees) to make a prediction. For more details concerning random forest, the reader is referred to the works of [START_REF] Breiman | Random forests[END_REF] and Friedman (2001).

Generalized additive models are generalized linear models where the response variable (Y i ) depends linearly on a sum of smooth functions of covariates [START_REF] Hastie | Generalized additive models[END_REF]. The general notation is given below:

g(µ i ) = A i θ + j 1 f j (x ji )
where g(.) is a link function,

µ i = E(Y i ) and Y i EF (µ i , φ) with EF (µ i , φ)
denoting an exponential family with mean µ i and scale parameter φ. A i is a row of the model matrix for any parametric model component (intercept) and θ the corresponding parameter vector. Finally, f j are smooth functions of the covariates x j . By using smooth functions to represent the dependence between covariates and response variable, this model has a great flexibility. However, choices must be made especially about the type of smooth functions used and their associated smoothness. An extended review on the smooth functions can be found in [START_REF] Wood | Generalized additive models: an introduction with R[END_REF]. The visualization of the smooth functions associated with each covariates is really informative and can bring insight on the response variable. In the same spirit of the previous section, the permutation feature importance can be computed with the vip package.

Fitting the models

The different models tested are presented in Table 2.3. The linear models are fitted with the lm function in base R, GAM models are fitted with the mgcv package and the random forest with the ranger package. To assess in a objective manner the performances of the formulas, GAMs and random forest we perform a repeated (10 repetitions) of a 10 fold cross-validation. The RMSE reported in the results section are the average of the out-of-sample RMSE.

Results

Comparison of the models

The averages of the out-of-sample RMSE are presented in the Figure 2.6. The model with the highest RMSE is the simplest model: the linear model with no transformation of the covariates. It is followed by the linear model with the log transformed variable, then the GAM. The lowest RMSE is obtained by the random forest. This observed order was expected as more complex algorithms are more able to learn complex patterns and therefore yield better results. The permutation feature importance is computed for the best algorithm (random forest) fitted on all data and is shown in Figure 2.7. It is clear that the significant wave height has the most predictive power for this model. It is followed by the wave period, the slope of the beach and the wave directional spreading. The least predictive variable is the wave direction. decrease. This observation is not logical as higher wave should lead to larger runup elevation on the beach. This could be due to the lack of observations during these conditions (high waves). The smoothing function of T p is monotonic and positive. This is expected as waves with larger period lead to increased runup on the beach. On the contrary, the smoothing function of the wave directional dispersion is monotonic negative meaning that waves with lower directional dispersion lead to larger value of Runup. Finally, the slope of the beach seems to have a slight positive influence on the Runup elevation and the wave direction does not seem to be a good predictor as it was excluded by the regularization.

The variable importance computed by the feature permutation method is shown in Figure 2.10. Such as the random forest, the most predictive variable is H s . The second most predictive variable is wave directional dispersion, it is

Predictions on buoy data

To get an idea of the predictions of our models (linear model with log transformation, GAM and random forest), we decided to predict the runup on the beach for each observation at the wave buoy. Because the slope of the beach is not available for each observation it has been removed from the models. The From this figure we can see that GAM does not seem to generalize well for unseen data. Indeed, some predictions of the runup are negative and the runup elevation seems to decrease as the significant wave height increases. The predictions seems more reasonable for the linear model or the random forest.

Discussion

In this section we will discuss the advantages and downsides of statistical learning methods in the modeling of the runup.

Statistical models are a lot less computationally demanding than process based model such as Xbeach for run modeling. In addition, the statistical models are also easier to train and calibrate than the process based models and they provide tailored prediction for one site. They can bring valuable insights by analyzing the predictive power of each covariates used in the model. By using the GAM, the influence of each covariate can be visualized and analyzed.

Unlike process based methods, the statistical methods do not rely on physical processes. Consequently, they need a large amount of data to take into account such processes and this constitutes their main limitations. In predictive usage, a SLM trained on a short dataset will encounter numerous events outside the distribution the training data and the associated predictions will be erroneous. This is the case for the GAM model where it predicts negative runup during unseen conditions.

One solution could be to collect more data, however, in the case of runup modeling on the GPB it does not appear as the best solution. Indeed, to improve the models presented in this work, data collected during energetic events (large H s and T p ) are needed. Such events could be dangerous for the operators measuring data. The cost of measurements campaigns can also be a hindrance.

A potential solution for the lack of data could be to switch to Bayesian framework in order to integrate prior knowledge of the physical process and to estimate uncertainties. When there are unseen conditions, the model will rely more on the prior knowledge than in data leading to more parsimonious predictions.

Conclusion

We compared the performances of linear models, random forest and generalized additive model to find the best predictive model for the runup elevation on the Grande Plage of Biarritz. This comparison was made with a repeated 10-fold cross validation to assess objectively the performance of these methods on unseen data despite the low number of observations.

The results and variable importance were analyzed for two models: random forest and GAM. Random forest yielded the best predictions with a RMSE of 19.3 cm. The variable importance analysis showed that the variables with the most predictive power were the significant wave height and wave period at peak. For the GAM model, the RMSE was slightly larger (20.5 cm). Even though it was not the best predictive model, GAM gave valuable insights with the visualization of the effect of each covariate on the prediction (Smooth functions). In addition,the importance variable analysis for this model was slightly different as the two most predictive variables were the significant wave height and the wave directional spreading.

For the wave runup prediction, SLM are a great alternatives to empirical formulas or numerical modeling. By learning directly from the observations they provide "tailored" predictions for one site. In addition, they can provide valuable insights by studying the predictive power of the covariates or by looking at the relationships (smooth functions) between the covariates and the response variable in case of generalized additive models.

Wave Climate Characterization with statistical downscaling

Waves play a major role in many coastal processes ranging from sediment transport process to coastal flooding process. The influence of the waves on sediment transport and erosion processes is discussed in the work of [START_REF] Munk | Refraction of ocean waves: a process linking underwater topography to beach erosion[END_REF] or [START_REF] Masselink | Introduction to coastal processes and geomorphology[END_REF]. Waves also have an impact on the coastal flooding process. This has been demonstrated in the first part of this chapter where wave characteristics showed a great predictive power for a storm surge model and also in the work of Bertin et al. (2015) where it has been proven that the waves increased significantly the storm surge level during Xynthia and Joachim storm events.

Deep knowledge about sea state characteristics and its seasonality is therefore essential not only for coastal management but also for safety at sea and on the coast. The wave climate can be characterized by analyzing wave data which can originate from two sources: observations (wave buoy, satellites) or modeled data (numerical models).

In the framework of the MICROPOLIT project, a preliminary work has been done to characterize the wave climate on the study site of Anglet [START_REF] Callens | Analyse des forçages physiques et climatiques influant sur le transport des micropolluants dans l'estuaire de l'Adour[END_REF] 1 . The wave data at the location of the buoy of Anglet come from a simulation covering the period from 1949 to 2014 with the spectral wave model WWII [START_REF] Roland | A fully coupled 3D wave-current interaction model on unstructured grids[END_REF]. Temporal decomposition was performed to extract the trend and seasonality of H s and extreme value analysis was performed to investigate the potential value of extreme H s and their return period. For H s a slight increase was observed over the years and a strong seasonality was highlighted with larger waves during the winter.

In this last part of chapter 1, we aim to characterize more precisely the wave climate at the Anglet's buoy location by studying the relationships between the wave characteristics observed at the buoy and the synoptic weather over the North Atlantic ocean. To reach this objective, we employ the statistical downscaling method which tries to find statistical relationship between a local predictand (wave characteristics observed at the buoy) and a regional scale predictor (SLP field) which has been categorized into different weather types to facilitate the analysis.

Data

Sea level pressure Sea level pressure data were provided by the research data archive of Colorado2 . They come from two reanalysis:

• CFSR (Climate Forecast System Reanalysis) from 1979 to 2010 with 6 hour time-step with 0.5 • resolution [START_REF] Saha | Ncep climate forecast system reanalysis (cfsr) 6-hourly products[END_REF].

• CFSv2 (Climate Forecast System Version 2) from 2010 to 2020 with 6 hour time-step with 1 • resolution [START_REF] Saha | Ncep climate forecast system version 2 (cfsv2) 6-hourly products[END_REF].

The spatial points of the first reanalysis have been sampled to correspond to the spatial points of the second reanalysis which has a coarser resolution.

The spatial domain of the data was restricted to North Atlantic: from 20 • N to 60 • N and from 70 • E to 20 • W with a spatial resolution of 1 • . In land points (red points in Figure 2.12) are not considered for this work as they show a greater variability in atmospheric pressure [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF].

Wave characteristics

The wave parameters such as significant wave height H s , peak period T p , di-

Statistical downscaling

The statistical downscaling method aims to find statistical relationship between local predictands (wave characteristics) and regional-scale predictors such as sea level pressure and sea level pressure gradient [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF].

This method relies on a weather pattern approach: it uses several data mining techniques such as dimensionality reduction and clustering to simplify the atmospheric pressure data into weather types, then it finds statistical relationship between wave characteristics and each weather type. This method allows for the study of seasonal and inter-annual variability of the sea state [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF]. The simplified methodology is as follows:

1. Processing of the SLP data

• Computation of the squared SLP gradient

• PCA to reduce dimensionality of the SLP and SLP gradient data

• Clustering into different weather types (WT)

2. Analysis of statistical relationships between weather type and sea state 3. Forecast of the sea state by using the statistical relationships found previously

Computation of the SLP gradient

The first step of this method is to prepare the regional scale predictors: sea level pressure (SLP) and the squared gradient of SLP (SLPG). The gradient of SLP is important for this analysis as it is a better representation of the winds responsible for the wave generating process [START_REF] Wang | Climate change signal and uncertainty in projections of ocean wave heights[END_REF]. SLP and SLPG are preferred over wind variables because global models are known to better represent sea level pressure than winds variables [START_REF] Caires | Projection and analysis of extreme wave climate[END_REF].

The squared gradient is computed for each point from the values of the four nearest grid points using the weights proportional to the inverse of the distance [START_REF] Espejo | Spectral ocean wave climate variability based on atmospheric circulation patterns[END_REF]. It is the sum of the squared zonal and squared meridional SLP gradients, which is proportional to wind energy [START_REF] Wang | Climate change signal and uncertainty in projections of ocean wave heights[END_REF].

Temporal averaging

The temporal range and the region of the predictor are really important when predicting local predictand. The generation time of the wave in the Atlantic

Ocean must be taken into account. For north Atlantic Ocean, it is less than 5 days [START_REF] Hegermiller | A multimodal wave spectrum-based approach for statistical downscaling of local wave climate[END_REF]. The ESTELA method [START_REF] Pérez | EStelA: A method for evaluating the source and travel time of the wave energy reaching a local area[END_REF] can be used to determine the number of days to take into account. This method allows one to know the region of wave generation and the temporal range needed for the predictor. In this work, we choose to average the SLP and SPLG data over 3 days (current day and 2 previous days) such as the work of [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF]; [START_REF] Espejo | Spectral ocean wave climate variability based on atmospheric circulation patterns[END_REF].

Dimensionality reduction

The SLP and SPLG data for each point averaged over 3 days for the period 1970-2020 are stored in a 15065 × 5026 matrix. The dimensionality of this dataset must be reduced for the next step which is the clustering with Self organizing maps (SOM). A PCA is applied on this data to reduce the dimensionality while conserving a large proportion of the variability of the data.

Only the n first components of the data representing 99% of the variability of the data is kept.

Clustering into weather regimes with SOM

The clustering with Self organizing maps (SOM) is done on the data projected on the n first components. SOM are a type of neural networks that rely on competitive learning. This algorithm aims to create a map of the data based on an iterative process involving neurons connected together by a grid. For one observation, the closest neuron is designated as the best match unit (BMU).

The coordinates of the BMU are then shifted toward the observation point at a designated rate (learning rate). The neighboring neurons of the BMU are also slightly shifted toward the observation. This process is repeated several times for all the observations of the training set. After the training process, the grid of neurons are approximating the data distribution and new observations can be assigned to a neuron (best matching unit) which can act as a clustering.

To cluster the projected data, other methods can be employed such as k-means algorithm (KMA) or maximum dissimilarity algorithm (MDA). A complete comparison of the three algorithm are presented in the work of [START_REF] Camus | Multivariate Wave Climate Using Self-Organizing Maps[END_REF]. In this work, SOM was chosen over the two other algorithms as it is the best technique to graphically characterize the multidimensional wave climate [START_REF] Camus | Multivariate Wave Climate Using Self-Organizing Maps[END_REF]. By projecting the classification on a two bi-dimensional lattice with proximity information (similar groups are closer), this method allows the visualization of patterns in multidimensional data which simplify the analysis of such data.

The clustering with SOM algorithm is performed with the R package kohonen.

Results

Weather types with SOM clustering PCA was performed to reduce the dimensionality of the data before the clustering. Among all the components, the first 26 were kept for the clustering

Relating weather types to sea state

The wave climate can be characterized by looking at the distribution of the wave parameters (H s , T p , θ) depending on the WT. winter in terms of wave characteristics. The WT 24 display a high-pressure system located on Europe while the WT 50 display a low pressure system located on England. There are significant differences in wave characteristics between these two WT. Indeed, the WT 50 (most frequent in winter season) display higher values and variability for significant wave heights. On the contrary, the WT 24 display higher values and variability for peak wave period.

There is also a small difference in wave direction between these two WT: the mean wave direction is northwest for WT 24 and west-northwest for WT 50.

Discussion

The work presented in this short article is mainly exploratory and was aimed to characterize the wave climate and its seasonality at the Anglet's buoy location.

However, the characterization of wave climate at one location is only one of the possibilities offered by the statistical downscaling method. Many extensions of this work can be considered as this method can be employed to:

• reconstitute the wave spectrum daily [START_REF] Espejo | Spectral ocean wave climate variability based on atmospheric circulation patterns[END_REF] or monthly [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF] • perform wave climate projections by using SLP fields from different climate change scenarios [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF][START_REF] Perez | Statistical multi-model climate projections of surface ocean waves in Europe[END_REF] • perform extreme value analysis for wave parameter depending on each weather type [START_REF] Camus | An atmospheric-to-marine synoptic classification for statistical downscaling marine climate[END_REF][START_REF] Rueda | A multivariate extreme wave and storm surge climate emulator based on weather patterns[END_REF] When we use this method, we make 3 important hypothesis [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF]:

• Variability of the local variable should be explained by the statistical connection

• Changes in the mean climate should lie within the range of its natural variability

• The relationships should be stationary.

A main limitation of statistical downscaling method lies in the fact that these hypothesis might not hold true in a context of global change [START_REF] Camus | A weathertype statistical downscaling framework for ocean wave climate[END_REF]. Indeed, extreme events that did not occur in the past may occur in the future and such events would lie outside the range of the past variability.

Moreover, global change may influence the statistical relationships observed in the past. It is worth noting that data availability is also a limitation for the statistical downscaling method. This method need a lot of data to see as much different events as possible. A site with low data availability can display large inhomogeneities when using this method [START_REF] Rueda | Marine climate variability based on weather patterns for a complicated island setting: The New Zealand case[END_REF].

Conclusion

The statistical downscaling method based on weather types allowed us to better understand the relationships between the synoptic weather in the North Atlantic Ocean and the wave characteristics observed at Anglet's buoy.

We showed that among the 100 weather types, some of them are representative of the mains circulation patterns in the North Atlantic Ocean such as the North Atlantic Oscillation (NAO) or the East Atlantic Oscillation (EA).

By studying the seasonal probabilities of occurrence of each WT, we highlighted large differences between seasons in terms of WT and variability. We also demonstrated that the distributions of H s and T p differ significantly between weather types. We observed that WT with low-pressure systems, more frequent in winter, are generally characterized by waves with larger H s and T p in comparison with WT with high-pressure systems or homogeneous SLP fields which are more frequent in summer.

For the disaster risk reduction (DRR) process, statistical downscaling method is an interesting tool for coastal researchers or stakeholders because it gives a better understanding of the wave characteristics and their variability at one location. As stated in the discussion section, the statistical downscaling method in this work is used only to characterize the wave climate. However, it can be extended to perform wave climate projections or extreme value analysis which can bring insights respectively on the evolution of the future risk and on climate variability of extreme events. 58

Conclusion

This chapter showed the interest of using SLM to model processes involved in coastal flooding. Indeed, in the first two applications of this chapter, we showed that they can predict new observations accurately when the new data is in the range of the training data (interpolation) and that they can provide knowledge about these processes or the study site with variable importance analysis.

In the last application of this chapter, we demonstrated that knowledge can also be acquired with unsupervised learning methods (PCA, SOM). These methods can be employed to detect patterns from big data and to gather similar observations into groups which is particularly useful in wave climate characterization.

Among the three examples of statistical modeling presented in this chapter, only the storm surge modeling can be considered for operational use in an EWS. Indeed, the statistical models trained to predict wave runup were not reliable for energetic conditions (storm events) due to a lack of data. Statistical downscaling could be used to predict wave characteristics depending on the synoptic weather type. However the wave characteristics can be more easily and accurately simulated with a spectral wave model. In this chapter, we show that SLM, especially ensemble methods (random forest and gradient boosting trees) can be employed to improve significantly the spectral wave model predictions at a specific location. In the literature, numerous methods combining statistical models and process-based models have already been proposed to improve the numerical wave model predictions (Babovic et al., 2001;Makarynskyy et al., 2005;Moeini et al., 2012;Deshmukh et al., 2016;Londhe et al., 2016). However, they all have in common to use shallow neural networks as off-the-shelves method to perform data assimilation. Unlike previous works on this subject, we compare the performances of neural networks with ensemble methods that have never been used for this task. In addition, a special attention is given to hyperparameter tuning. Hyperparameters are parameters associated to machine learning algorithms that must be chosen before the training. These parameters influence greatly the final performance as they control the structure or the learning of the algorithm.

In the literature, not a lot of importance is given to hyperparameter tuning and usually only a few combinations of hyperparameters are tested with grid search method. In this work, bayesian optimization, an automatic method, is employed to find the optimal hyperparameters. Finally, the potential benefits of data assimilation are investigated in a real case scenario: the computation of wave run-up on the Grande Plage of Biarritz.

Scientific output:

This work resulted in the publication of an article (presented below) and two communications.

Publication:

Abstract

The main objective is to present alternative algorithms to neural networks when improving sea state forecast by numerical models considering main spectral bulk parameters at a specific location, namely significant wave height, peak wave period and peak wave direction. The two alternatives are random forest and gradient boosting trees. To our knowledge, they have never been used for error prediction method. Therefore, their performances are compared with the performances of the usual choice in the literature: neural networks. We showed that the RMSE of the variables updated with gradient boosting trees and random forest are respectively 20 and 10% lower than the RMSE obtained with neural networks. A secondary objective is to show how to tune the hyperparameter values of machine learning algorithms with Bayesian Optimization. This step is essential when using machine learning algorithms and can improve the results significantly. Indeed, after a fine hyperparameter tuning with Bayesian optimization, gradient boosting trees yielded RMSE values in average 8% to 11% lower for the correction of significant wave height and peak wave period. Lastly, the potential benefits of such corrections in real life application are investigated by computing the extreme wave run-up (R 2% ) at the study site (Biarritz, France) using the data corrected by the different algorithms. Here again, the corrections made by random forest and gradient boosting trees provide better results than the corrections made by neural networks.

Introduction

Nowadays, numerical wave models are routinely used to forecast wind generated waves. Although they provide satisfactory predictions at a regional scale and during mean wave conditions, it has been shown that they are less accurate for forecasting at a specific location (Londhe et al., 2016) and have a tendency to underestimate wave height during energetic wave conditions. This underestimation has been observed in different state-of-the-art wave models: see the work of Arnoux et al. (2018) for Wind Wave Model II (WWMII) and WAVEWATCH-III (WW3) ; Rakha et al. (2007) for WAve Model (WAM) and Moeini et al. (2012) for the Simulating WAves Nearshore model (SWAN). The errors in wave predictions are mainly due to inaccuracies in the wind input that forces the model. The winds used as forcing are numerically simulated and are known to underestimate high wind speeds (Moeini et al., 2012). This results in the underestimation of wave parameters by numerical wave models. Simplifying assumptions, approximations employed in the modeling process, discretization of the domain and a potentially wrong parametrization of the model can also be sources of inaccuracies in wave model predictions (Babovic et al., 2001(Babovic et al., , 2005)).

When observation data are available, data assimilation can be used to improve the predictions made by numerical models. There are 4 main categories of data assimilation procedures [START_REF] Refsgaard | Validation and intercomparison of different updating procedures for realtime forecasting[END_REF]Babovic et al., 2001): updating the input parameters, updating the state variables, updating the model parameters and finally updating the output parameters. The last procedure is called "Error prediction" method and is the most suitable approach to improve model predictions of different output variables at a specific location (Babovic et al., 2005). This procedure presents several advantages comparing to the other data assimilation procedures. First, it covers inaccuracies coming from all sources because it improves directly output variables. In addition, it can use a combination of external variables such as meteorological or wind data to increase the accuracy of the predictions. Lastly, it is easy to implement because it consists in only three steps and does not require multiple runs of numerical wave model. First, the deviations between the modeled values and measured values are computed. Then, machine learning algorithms are used to forecast these deviations. Finally the deviations predicted by the algorithms are incorporated to the predictions of the numerical model for the next time steps, resulting in a more accurate wave forecast.

This method has been successfully applied on hindcast data (Makarynskyy * Corresponding author et al., 2005;Deshmukh et al., 2016) and has even been implemented in real time setting in the works of Babovic et al. (2001) and Londhe et al. (2016).

To our knowledge, only artificial neural networks have been tested to forecast the errors in the data assimilation. However, according to the so-called "No Free Lunch" theorem, there is no single model that works best for all problems [START_REF] Wolpert | The supervised learning no-free-lunch theorems[END_REF]. It is therefore necessary to try multiple models and find the one that works best for our particular problem. The performance of artificial neural networks must be compared with other algorithms in the data assimilation task. Random forest and gradient boosting trees are strong candidates for this comparison. Indeed, these two methods are known for their performance and unlike neural networks, they also provide valuable information by computing the predictive power of each variable used as input. The predictive power or variable importance refers to how much a model relies on that variable to make accurate predictions. A variable with high predictive power means that its values have a significant impact on the prediction values. By contrast, a variable with low predictive power have a limited impact on the prediction values and it can be substracted from the model to make it simpler and faster.

To explore the performance of random forest and gradient boosting trees, we use as a test case the Basque coast (South west of France). Every winter, the basque coast faces numerous coastal flooding events. To prevent and mitigate the risk of flooding, wave forecast are used to compute the extreme run-up values either by using parametric models such as the formula of Stockdon et al. (2006) or process based models such as Xbeach (Vousdoukas et al., 2012;de Santiago et al., 2017). In both cases, the accuracy of this forecast is of utmost importance as the issuing of the early warning depends on it, especially during energetic wave climate where coastal flooding risk is the highest. In this study, we employ the error prediction method with the different machine learning algorithms and use local meteorological conditions and measured wave parameters from a local buoy to improve the wave forecast. Lastly, we investigate the potential benefits of using such corrections in the computation of extreme run-up values.

This study aims to present two alternatives (random forest and gradient boosting trees) to neural networks by comparing their performances when improving regional numerical models. A secondary objective is to show how to tune the hyperparameter values of machine learning algorithms with Bayesian Optimization. In machine learning, a hyperparameter is a parameter whose value is specified by the user before the learning process begins, it will affect how well a model trains and therefore it will have a non negligible impact on the final results. Bayesian optimization is an efficient hyperparameter optimization algo-rithm and it is widely used to optimize the results of any given machine learning method.

Lastly, we investigate if the error prediction method makes a difference in a real application such as the computation of extreme run-up for the beach of Biarritz. Section 2 will introduce the study area, the data and all the statistical methods used. Results will be presented and discussed in Section 3. Finally, Section 4 will cover the conclusion.

Data and Methods

Study site and Data

The Basque coast is a 150 km long rocky coast facing the Bay of Biscay (Figure 1). Every winter, it is battered by numerous storm events. This results in frequent and sometimes intense coastal flooding which can severely damage seafront infrastructures. The city of Biarritz is particularly affected as the buildings and infrastructures are located right behind a sea wall that is located at the top of the beach. The damages associated with coastal flooding are costly for nearshore cities which try to prevent and mitigate the risks by developing early warning systems. Such systems rely on the knowledge of the sea state and its forecast. This work focuses on the forecast improvement of three wave integrated parameters which describe the sea state in this area: the significant wave height (H S ), the peak period (T p ) and the peak wave direction (θ p ). Direct measurements of these parameters are obtained from the National Center for Archiving Swell Measurements (L'her et al., 1999). They were made by a directional wave rider buoy (DWR MKIII) operated by the Centre for Studies and Expertise on Risks, Environment, Mobility, and Urban and Country Planning (CEREMA) and the University of Pau and Pays de l'Adour (UPPA). The buoy is located a few miles off the Basque Coast (Figure 1) at 50 meters water depth. Since its deployment in 2009, this buoy have been recording the parameters of interest every 30 minutes. The measuring range of this buoy is [-20m; 20m] for heave motion, [1.6s; 30s] for wave period and [0 • ; 360 • ] for wave direction. It has a resolution of 1 cm in heave motion and a directional resolution of 1.5 • . To be consistent with the numerical wave data and meteorological data, a 1 hour time step was adopted for the buoy data.

The three parameters simulated at the buoy coordinates by the Meteo-France WAM model were provided by the Copernicus Marine Environment Monitoring Service. This reanalysis ("ibi reanalysis wav 005 006") covers the period 2007-2019 with a hourly time-step. The MFWAM model is derived from the third generation wave model WAM [START_REF] Group | The WAM model-A third generation ocean wave prediction model[END_REF]. It is forced by wind fields obtained from a regional numerical weather prediction model (AROME). A more complete description of the MFWAM model can be found in Lefèvre and Aouf (2012).

Meteorological data, including average wind speed above 10 meters, wind direction and atmospheric pressure were furnished by the French national meteorological service MeteoFrance. The data were collected hourly by the meteorological station of the Biarritz airport, located only a few kilometers from the study site (Figure 1). It covers the period ranging from 2013-01-01 to 2018-12-31. By assembling the wave buoy data, the wind wave parameters and the meteorological data we obtain a dataset of 41439 hourly observations ranging from 2013-01-01 to 2018-12-31. In this work, we are improving the wave forecast by correcting the systematic errors of the wind wave model. Therefore, we are not considering any temporal effects while improving H S , T p and θ p . The dataset was randomly divided into 2 parts: the training part containing 70% of the observations (n = 28797) and the testing part containing the remaining 30% (n = 12342).

Error prediction method

The error prediction method consists in three steps:

• Step 1: Deviations between model predictions and measured values are computed:

E model = X measured -X modeled ,
where E model is the error of the model, X measured is the measured value of an output variable provided by the wave buoy and X modeled is the value of the same variable computed by the wave model.

•

Step 2: E model is predicted with an appropriate supervised machine learning algorithm.

• Step 3: The predicted error is added to the prediction of the wave model to obtain an updated numerical prediction:

X updated = X modeled + E predicted ,
where X updated is the updated prediction of wave model and E predicted is the predicted error given by the supervised learning method.

This method is repeated separately for each output variable to improve (H s , T p , θ p ). The performance of this data assimilation method relies on two things: the quantity of data and the machine learning algorithm used. Since the machine learning algorithm are generally more suited to interpolate rather than extrapolate, the available data for learning process should cover as much as possible the range of all the probable events in the study area. Concerning the learning method, only neural networks have been used for the step (2) of the error prediction method to our knowledge (Makarynskyy et al., 2005;Moeini et al., 2012;Londhe et al., 2016). Because we want to compare the performance between different machine learning algorithms, we use random forest and gradient boosting trees. All the tested algorithms use the same input variables to improve the model accuracy: the three wave parameters (H S , T p , θ p ) given by the numerical model, the atmospheric pressure, the wind direction and speed.

Neural networks

Artificial neural networks have been extensively used in the domain of wave modelling (Deo et al., 2001;[START_REF] Makarynskyy | Artificial neural networks in the forecasting of wave parameters[END_REF]Makarynskyy, 2005;Mandal and Prabaharan, 2006) or wave parameters assimilation (Makarynskyy et al., 2005;Moeini et al., 2012;Londhe et al., 2016). It is why technical details will be avoided in this study and only the general concepts will be presented. The readers can find more details and information on the working of neural networks in [START_REF] Liang | Neural network fundamentals with graphs, algorithms and applications[END_REF] or [START_REF] Friedman | The Elements of Statistical Learning[END_REF].

The most common class of neural networks is the multilayer perceptron. The neurons in this network are organized in three layers: the input layer that receive the input variables, the output layer that performs the final predictions and between these two layers there is the hidden layer. Neurons in the hidden layer transmit the signal to the output layer by transforming the weighted sum of the neurons present in the input layer with a non linear function called activation function. The weights between each neuron of the network are adjusted through the iterative process of backpropagation to minimize the error between the variable we want to predict and the variable predicted by the network (output layer).

As other machine learning methods, hyperparameters need to be specified before the training of neural networks. Some hyperparameters control the network architecture (number of neurons, layers, activation function used, etc...) while others control the training process (learning rate, bacth size, number of epochs, etc...). Hyperparameters must be tuned carefully in order to achieve optimal results with neural networks.

Tree based algorithms

Unlike neural networks, random forest and gradient boosting have never been used in the error prediction method. They are state-of-the-art ensemble learning techniques for classification and regression tasks. An ensemble learning technique commonly refers to a method that combines the predictions from multiple machine learning algorithms, called base learners, to produce more accurate predictions.

Random forest is an algorithm that builds many decision trees in parallel. These trees are the base learners for random forest and they have the following characteristics:

• Each tree is built using a different bootstrap sample of the data-set. This mechanism is called bagging.

• At each node, a given number (hereafter "mtry") of variables are randomly sampled as candidates at each split. The best split point is then selected within this random set of variables. This process is called feature sampling.

The value "mtry" is fixed before growing the forest.

• Unlike the classification and regression trees of [START_REF] Breiman | Classification and regression trees Chapman & Hall[END_REF], the trees in random forest are fully grown (no pruning step).

Bagging and feature sampling are the core principles of random forest. They are two randomizing mechanisms which ensure that the trees are independent and are less correlated with each other. The final prediction of a random forest is obtained by averaging the results of all the independent trees in case of regression or using the majority rule in case of classification.

The most important hyperparameters in random forest are the number of trees and "mtry": the number of variables randomly sampled as candidates at each split when building the trees.

Gradient boosting is an algorithm that trains many weak learners sequentially to provide a more accurate estimate of the response variable. A weak learner is a machine learning model that perform slightly better than chance. In case of gradient boosting trees, the weak learners are shallow decision trees. Each new tree added to the ensemble model (combination of all the previous trees) minimizes the loss function associated with the ensemble model. The loss function depends on the type of the task performed and can be chosen by the user. For regression, the standard choice is the squared loss. By adding sequentially trees that minimize the loss function (i.e. follow the gradient of the overall loss function), the overall prediction error decreases. Technical details about gradient boosting trees can be found in (Friedman, 2001).

Many hyperparameters have to be tuned for gradient boosting trees, some of them control the gradient boosting process, such as the learning rate, the number of trees to be used whereas others regulate the construction process of the trees: minimal node size, sample of the dataset to be used, maximum depth.

Hyperparameter tuning

Hyperparameters influence significantly the training of the machine learning algorithms and therefore the quality of their predictions. The objective of hyperparameter tuning is to find the values of hyperparameters that yield the lowest error (RMSE in our case) for unseen data. Two types of methods exist to find the optimal values of hyperparameters: uninformed or informed.

In uninformed methods, many combinations of hyperparameter values are tested one after the other and the best combination is the one that yields the lowest error on unseen data. The values of hyperparameters are either sampled randomly (random search) or sampled along a grid (grid search). In both cases, each combination tested are independent from another. With grid and random search, it is not guaranteed to find the optimal set of hyperparameters and it usually requires a lot of iterations (combinations tested).

In informed methods, the results obtained by the past combinations are used to choose the next combination to evaluate. Bayesian optimization algorithm is an informed method that aims to minimize an objective function, in our case the errors of the machine learning algorithms on unseen data. First, it builds a probability model (Gaussian process) of the objective function. Then it uses this surrogate model to select the most promising values of hyperparameters to evaluate. Once the promising combination of values have been evaluated, the probability model is updated and searched again for the most promising combination. This process is repeated several times. This method is employed in this article because it is very efficient for tuning hyperparameter values and it usually requires less iterations than uninformed methods (Bergstra et al., 2011). In-depth details of this method are given in the works of Snoek et al. (2012); Marchant and Ramos (2012) and Shahriari et al. (2015).

Training the algorithms

The machine learning algorithms described above are trained to predict the deviations of H s , T p or θ p (one model for each variable), using 6 input variables: the three wave parameters (H S , T p , θ p ) given by the numerical model, the atmospheric pressure, the wind direction and speed.

The neural networks are built and trained with the R package keras. The input variables are centered and scaled to improve the result of neural networks and the weights are updated with the adam optimization algorithm [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. Random forest and gradient boosting model are fitted in R using respectively the ranger package which provide fast implementation of Random Forests (suited for high dimensional data) and the xgboost package which is an efficient R implementation of the gradient boosting framework from [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF]. The input variables are not centered or scaled before the training of random forest and gradient forest because it does not influence the training of these algorithms.

The training is done twice: once with the default values of the hyperparameters in the R packages and once with the optimal values found with the Bayesian optimization method.

The best hyperparameter values are found by the means of Bayesian optimization method coupled with a 5-fold cross validation in the training dataset. That is, the training data are split into five equal-sized partitions and a machine learning model is recursively built on four partitions (80% of the training data) with a given hyperparameter combination. A performance metric, in our case the root mean square error is assessed on the remaining partition (20% of the training data). The resulting five performance metrics are averaged to provide an estimated out-of-sample performance of the respective hyperparameter combination. The objective function to minimize for the Bayesian optimization method is the average out-of-sample performance value. The Bayesian optimization for our data is performed using the R package RBayesianOptimization. First, random combinations of hyperparameter values are evaluated to serve as search base for the informed method (5 in this study), then an acquisition function (upper confidence bound) is used to find the next combination values to evaluate (this step is repeated 25 times).

Results and Discussion

Model comparison

To assess the accuracy of the numerical model and the proposed corrections, several metrics are computed including the root mean square error (RMSE), the correlation coefficient, the bias and the scatter index (SI). The bias represents the average error between the observed and modeled data and allows one to detect under or over estimation of the value of one parameter. The scatter index is a measure of the error normalized by the observation values. It is a standard metric for wave model inter-comparison (Londhe et al., 2016). More details about the computation of these two metrics can be found in the work of [START_REF] Mentaschi | Problems in RMSE-based wave model validations[END_REF]. The metrics are computed twice: once with the whole test data and once with a subset of the test data where H s > 3m because the underestimation of H s is known to become larger above this height (Arnoux et al., 2018).

Table 1 presents the metrics obtained with no assimilation (numerical model) and after a preliminary data assimilation with the three machine learning algorithms. The term "preliminary" refers to the lack of hyperparameter tuning. The learning has been done using the default hyperparameter values given in table 2.

For significant wave height, the numerical model shows a negative bias. This indicates that the MFWAM model has a tendency to underestimate H s such as other wind wave models (Moeini et al., 2012;Arnoux et al., 2018). The negative bias increases as the value of H s becomes larger (H s > 3m), meaning that H s is more likely to be underestimated during energetic events. For the peak period, the numerical model shows a positive bias. When H s > 3m, the bias and the RMSE for this parameter are smaller. The predictions of T p are therefore better during energetic conditions. For wave direction, a small bias is observed in average and is greater when the waves are larger. The large difference in RMSE computed with data where H s > 3m and with all data is explained by the distribution of the wave direction according to the wave height. When the significant wave height is below 2 meters, the wave direction at the buoy is more variable (Figure S1, supplementary material) and the spectral wave model has more difficulties to predict correctly the direction. This can be confirmed by looking at the θ p errors of the numerical model: we see that they are larger and occur more often when H s < 2m (Figure S2, supplementary material). A potential explanation of this phenomenon could be that below 2 meters, the sea state is more likely to be influenced by local wind conditions which are difficult to reproduce by the spectral wave model [START_REF] Rascle | A global wave parameter database for geophysical applications. part 2: Model validation with improved source term parameterization[END_REF]. When the significant wave height is above 2 meters, the wave directions are a lot less variable and the predictions of the spectral wave model are more accurate.

When we look at the metrics computed with all data, we see that the correction made by the three machine learning algorithms removes the bias and greatly reduces the RMSE and the scatter index for H s and T p . For θ p , the mean bias is slightly larger after data assimilation for all algorithms. The correction of the machine learning algorithm could be less efficient for θ p due to the high variability of the observed deviations they try to model (see the explanation in the paragraph above). However, lower value of RMSE and scatter index and larger correlation coefficients still indicate that the corrected data are closer to the observed values at the buoy.

For the metrics computed with data where H s > 3m, the correction does not remove the bias for H s and T p but reduces it greatly. For the wave direction, the updated parameters are closer to the reality. Indeed, bias and RMSE obtained by the corrections are smaller than the numerical model and the correlation coefficients are larger for corrected data.

For this preliminary assimilation, random forest yields the best results for all the parameters. It reduces the RMSE values computed with all test data by 37.7%, 30% and 29% respectively for H s , T p and θ p . Gradient boosting trees is Note: h corresponds to the number of input variables (6 in our case).

close second and decreases the RMSE values by 33%, 24.5% and 22.8%. Finally, data assimilation with neural networks decreases the RMSE of H s , T p and θ p by 23%, 12.8% and 10.2%. As stated earlier, the performance of machine learning algorithms might depend on the choice of the hyperparameter values. The Bayesian optimization was therefore performed and optimal values were selected (Table 2). The selected hyperparameter values are quite different from the default values. Indeed, for neural networks, the best results were obtained with more epochs and more neurons in the hidden layer. For random forest, only the number of trees seems to have some effect on the results and models with a large number of trees performs better. Finally, for gradient boosting trees, models with a large number of trees and a small learning rate are preferred.

Metrics calculated with data corrected by the tuned machine learning algorithms are presented in Table 3. Overall, tuning the hyperparameter values has improved the results of all the algorithms. However, the degree of improvement differs depending on the algorithm. We observe the smallest improvements for random forest where the RMSE of every parameters seems to decrease by less than 1% in average. For neural networks, tuning hyperparameter values has a more significant effect by reducing the RMSE by 2 to 3% in average. The largest effect of tuning the hyperparameters are observed with gradient boosting trees. The RMSE is 8 to 11% lower for every parameter. The only exception is θ p computed with all data where we have a small increase (2%) of RMSE. In general, the mean bias for H s , T p and θ p remains the same before and after hyperparameter tuning expect for the bias of H s computed when H s > 3m which is significantly For this dataset, gradient boosting algorithm shows the best performances for all parameters. Assimilation with this algorithm decreases the RMSE values computed with all test data by 39.8% for H s , 33% for T p and 31% for θ p . For H s and θ p , the reduction are even lower for the RMSE values computed with H s > 3m: 47% for the significant wave height and 40% for wave direction. The performances of random forest for H s , T p and θ p are slightly better than the results obtained before tuning the hyperparameters: respectively 38.3%, 30.9%, 30.1%. The performances are also better for neural networks after hyperparameter tuning: it decreases the RMSE values by 24.8% for H s , 15.5% for T p and 12.5% for θ p . The differences in efficiency between neural networks and ensemble learning techniques could be explained by the architecture chosen for the neural networks. Indeed, this work shows the results for multilayer perceptrons with only one hidden layer which is the typical choice in the literature (Londhe et al., 2016;Moeini et al., 2012). By choosing an architecture with more hidden layers, the networks might be able to model more complex phenomena and bring a better improvement for the three wave parameters.

The distribution of the errors after the different corrections are presented in the figure 2. For all wave parameters, the distributions of the errors after a correction have narrowed and are now more centered in zero. The differences in performance between algorithms are confirmed with these violin plots. Indeed, when the correction is made with random forest or gradient boosting trees, the distributions of the errors are more narrow than the distributions of the errors obtained with neural networks. The difference in efficiency between random forest and gradient boosting trees is not distinguishable graphically. It is expected as the metrics of the two algorithms only differ by a few percents. For H s and T p , the corrections have also removed the bias observed for numerical model. The large errors of θ p for the numerical model (Figure 2) are observed when H s < 2m and are not corrected by the machine learning algorithms. Figures showing the observed values versus the corrected values are available in the supplementary material for the three wave parameters.

Predictive power of the input variables

In addition to their performance, random forest and gradient boosting algorithms can provide a measure of importance for each variable used as input. This importance indicates the predictive power of the variable. It can be used to sort variable from most to least predictive, allowing one to have more insight on the problem and to perform feature selection when there are too many input variables. The figure 3 shows the importance measure of each variable computed by the random forest depending on the parameter to improve. For H s and T p , the most important variables are the value of H s and T p modeled by the wind wave model. It is different for the direction where the most important variables are the value of θ p and H s given by the model. The predictive power of local meteorological variables is quite low, suggesting that local and instantaneous meteorological variables does not bring valuable information in the assimilation process. The wind wave formation process is not instantaneous and occurs in large regional scale, therefore using meteorological variables from the past (several days before) and from different locations (located in the ocean) could lead to a better predictive power which means better updated wave predictions.

Example of application

To investigate the potential effect of the different corrections in a real case scenario, the extreme wave run-up R 2% at the Grande Plage de Biarritz has been computed for the test period with the Stockdon formula (Stockdon et al., 2006) which uses H s and T p and the beach slope as parameters. The beach slope is fixed to 8% according to the work of Morichon et al. (2018). Using the extreme wave run-up calculated with the buoy data as reference, the metrics presented previously have been computed for the numerical model and the different corrections (Table 4). From this table, it is evident that the data corrected with machine learning algorithms provide wave run up values that are closer to the "real" values with lower RMSE, Scatter index and greater correlation coefficient. Although the bias remains, the correction made by the gradient boosting

Conclusion

In this work, random forest and gradient boosting trees were employed for the first time in the error prediction method. These ensemble learning techniques based on decision trees performed better than neural networks for improving the wave forecast of the Basque Coast. The correction made by gradient boosting trees yielded the best results for all the wave parameters: it reduced the RMSE values by nearly 40% for H s , 33% for T p and 31% for θ p . The reduction of RMSE values for random forest was only a few percents lower than gradient boosting trees. The corrections made by neural networks were significant but yielded reductions in RMSE not as high as the two ensemble learning techniques: 24.8% for H s , 15.5% for T p and 12.5% for θ p .

As expected, tuning the hyperparameters of the machine learning algorithms had a positive effect on the final results. However, the effect of the tuning differed depending on the algorithms. Indeed, random forest was less affected as it only reduced the RMSE values by 1% in average. The tuning had more effect on neural networks reducing the RMSE values by 2 to 3%. Gradient boosting tree algorithm was the most affected by hyperparameter tuning as the results were improved by 8 to 11% in average. One of the main advantage of random forest over gradient boosting trees is that it doesn't need this tuning step in order to yield great results. This is not negligible as hyperparameter tuning step can be time consuming and computationally demanding depending on the complexity of the search (number of hyperparameters).

Contrary to neural networks, Random forest and Gradient boosting trees provided valuable insights by giving the predictive power of each input variable. The predictive power of variable brings interpretability to the model and can give a better understanding of the variable we try to predict. For example, we know that the significant wave height modelled by the numerical wave model was the most important variable in the correction of the three parameters. In cases where there are a lot of input variables, knowing their associated predictive power helps developing more parsimonious models by keeping the pertinent variables and subtracting the less informative ones from the model.

The error prediction method has proven to be useful in improving wave forecast. This had an impact in a real life application by improving the accuracy of the extreme run-up computed at the Grande Plage de Biarritz. Here again the corrections brought by random forest and gradient boosting tree were better than the correction made by neural networks. The decrease in RMSE values was around 22% for the two ensemble techniques and 6.2% for the neural networks. Even though the differences in performance might not appear significant, it can make a difference when using these corrections in an early warning system. It is especially true when dealing with storm events where H s and T p are large.

The differences between machine learning algorithms observed in this article are specific to Biarritz site. The results might differ for another study site. Therefore, we can only advise to test and compare several machine learning algorithms to find the optimal one associated with the site of interest.

Finally, the assimilation made in this study did not account for the temporal aspect in the errors of the numerical model, it only corrected systematic errors of the wave model. In the future, this work could be extended by adding input variables containing temporal aspect. This could be the values of a modeled parameter at previous time steps such as the work of Londhe et al. (2016). In this framework, neural networks could perform better as they are known to handle efficiently time series. Other input variables could be also used to improve the wave forecast such as the meteorological data from the past or at different locations. Because the success of the error prediction method depends on the quantity of data, it would be also interesting to perform a sensitivity analysis on the quantity of data used in the training process. This could give us some insights on the minimal quantity of data required to obtain a desirable assimilation procedure.

Supplementary material: Using Random forest and Gradient boosting trees to improve wave forecast at a specific location

Observed data vs data corrected with different machine learning algorithms

Conclusion

In this chapter, we showed that SLM can be employed to improve significantly the predictions of numerical wave models at a specific location. We also proved the importance of comparing different SLM and choosing the optimal hyperparameters in order to obtain the best results. With at least several years of wave buoy data, this method can be easily implemented to improve the predictions of numerical wave models in an EWS.

The methodology presented in this chapter only corrects the systematic errors of the numerical model at a specific location. This work could be therefore extended by including temporality in the error prediction method. To that end, errors of previous time steps could be included in the explanatory variable or methods specialized in time series prediction such as recurrent neural networks (RNN) could be employed in the same spirit as [START_REF] Zhang | Ocean wave height series prediction with numerical long short-term memory[END_REF]. In this work, they demonstrate that RNN can be used to improve the predictions of a numerical wind wave model by integrating both the local data and the temporality in the errors of the numerical wind wave model. Their proposed methodology with RNN led to a better accuracy than the error prediction method performed with usual machine learning algorithms. The work presented in this chapter could also be extended by correcting the whole wave field (multiple locations) predicted by the numerical model. The objective of this chapter is to demonstrate that SLM can be employed to create automatically a storm impact database with images from video monitoring networks. Data about the impact of a hazard is usually the limiting factor when building a data-driven model translating tide, wave and meteorological conditions into impact on the shore. Unlike data about the tide, wave or meteorological conditions, data about storm impact are rare, sparse and mostly come from archives or insurances data. So far, no methods exist to collect routinely storm impact data. We propose to use Convolutional neural networks (CNN), which are deep learning methods, to classify the video monitoring images into three storm impact regimes which are categories of coastal flooding risk. Several CNN architectures and methods to deal with class imbalance are tested on two sites (Biarritz and Zarautz) to find the best practices for this classification task. Transfer learning is also investigated to facilitate the transferability of this method to new sites. Once trained, the CNN can predict the storm impact regimes of newly created timestacks, generating an incremental storm impact database.

. Introduction

Databases containing information on past storm characteristics and their impacts on the coast are essential for the disaster-risk-reduction process. They enable scientists and coastal stakeholders to better understand the storm hazard in a specific area, to identify potential trends, and most importantly to assess coastal risks (present or future) through their use in the development and validation of early warning systems [1,2].

In these databases, storm impact is mostly represented as a qualitative variable with different categories. The different categories of storm impact are called "regimes" and are defined according to the Sallenger's scale [3]. This scale was originally derived to classify storm impact intensity based on the relation between wave-induced maximum water level and topographic elevations of the different sections of a natural beach. Recently, this approach has been extended to the estimation of storm impact intensity at an engineered beach backed by a seawall [4].

Due to the extreme and episodic nature of storms, databases covering a long period of time are necessary. Observed data about storm characteristics such as tide, wave, and Due to the extreme and episodic nature of storms, databases covering a long period of time are necessary. Observed data about storm characteristics such as tide, wave, and wind are abundant and have been collected routinely for decades. In addition, numerous reanalyses and hindcasts are available for these variables. On the contrary, data on storm impacts are more sparse and mostly come from archives [5][6][7] or insurance data [8,9]. A few examples of storm impact databases are: the RISC-KIT database, which contains storm impact information for nine study sites in Europe [10]; the SurgeWatch database [5] for the UK; and a database for the Basque coast [7]. Even though archives and insurance provide information, there are some limitations including the heterogeneity, the incompleteness of the data sources, and the consequent amount of work needed [7]. A solution to routinely create a storm impact database could be to use images provided by coastal monitoring stations that are now widely used worldwide to survey and study coastal processes.

In recent decades, video monitoring systems have proven to be valuable assets in the study of the coastal zone due to their cost-efficiency and their ability to provide a continuous stream of data including intense storm conditions. Video monitoring systems are generally composed of one or several cameras operated by a monitoring software such as Argus [11], HORUS, Kosta (www.kostasystem.com, accessed on may 2021) or Sirena [12]. The reader is referred to the work of Nieto et al. [12] for a comparison of the cited monitoring systems. These systems generate different types of images that can be applied to study coastal processes such as beach morphology changes, wave runup, and coastal currents [13,14]. Among the different types of images generated by the video monitoring system, timestacks images represent the time-varying pixel intensities along a particular cross-shore transect in the camera's field of view. They are used to perform wave runup parametrization [15][16][17], wave breaking detection [18], or intertidal topography [19] and also to estimate wave characteristics [20,21], sea level [22], and bathymetry [23,24]. Timestacks images have also been employed in the study of storm impact. In the work of Thuan et al. [25], they quantify the impact of two typhoons on the longshore-averaged shoreline changes based on the analysis of a series of timestacks images. To our knowledge, timestacks have not been used to directly measure storm impact regimes as defined previously. Image processing techniques are usually employed to transform the information contained in the images into quantitative measurements (runup elevation, wave height, shoreline). In this article, we propose extracting storm impact regimes (qualitative data) directly from the timestacks.

The storm impact regimes can be extracted from timestacks using two methodologies. The first methodology can be qualified as deterministic; it relies on image processing techniques and consists of two steps. First, the water line position is found by segmenting the image. The storm impact regime is then deduced by comparing the position of the waterline with the position of the defense infrastructure in the timestack. Different methods can be used to extract the waterline from timestacks images [26,27]. For example, Otsu's method [26] divides the pixels into two groups depending on their intensity values. It is not always robust and depends on the quality and lighting of the images. Most of the time, it requires rigorous and tedious human verification and correction [16]. The second methodology, presented in this article, relies on deep learning with convolutional neural networks (CNNs). CNNs are a class of deep neural networks, specializing in imagery analysis, that perform well on specific problems such as image classification and segmentation. First, timestack images are classified into storm impact regimes by human operators. Then, the CNN is trained to classify timestacks into storm impact regimes using the annotated dataset. During the training process, the CNN learns to simultaneously classify the images and which features to detect in order to achieve the best classification accuracy. Once the neural network has learned on the training dataset, it can be used to routinely analyze the timestacks produced by the video monitoring system and therefore create incrementally a qualitative storm impact database. This second methodology, based on a self-learning algorithm (CNN), allows for more automation compared to the first methodology because it does not require site-specific calibration [17].

has already been employed in coastal engineering domain [14,[START_REF] Valentini | Assessment of a smartphone-based camera system for coastal image segmentation and sargassum monitoring[END_REF] and usually results in faster training and better accuracy. In the case of storm impact recognition, where images of extreme storm impact regimes are rare by nature, this method can significantly improve the performances of CNNs. Moreover, it is reasonable to think that knowledge acquired at one site can be used to improve the performances on another site. This could be a non-negligible asset for the application of the method to a new site.

This paper aims to demonstrate the high potential of CNN methods to constitute a storm impact database using timestacks images provided by coastal video monitoring stations. Different methods are tested using images collected at two study sites. The best practice and the transferability of knowledge gained at one site to another are studied. In the following sections, the study sites and the video dataset are first described. The main features of the CNN implementation procedure are then shown in Section 3. Results and transferability of the CNN between the study sites are presented in Section 4 and discussed in Section 5. The main results are finally presented in the conclusion, Section 6.

Study Sites and Data

In this study, the storm impact intensity is classified into three storm impact regimes (Figure 1) derived from the Sallenger's scale [3]. The following three categories have been adapted for the timestack images:

• Swash regime: all the waves in the timestack are confined to the beach; • Collision regime: at least one wave in the timestack collides with the bottom of the seawall ; • Overwash regime: at least one wave in the timestack completely overtops the seawall. The CNNs were trained on timestacks images collected by video monitoring stations operating on two sites along the basque coast (Figure 2), namely the Grande Plage of Biarritz (GPB), and the Zarautz beach (ZB). The use of two sets of data acquired from sites with different geological and morphological characteristics and distinct responses to oceanic forcing makes it possible to assess objectively the ability of CNN to detect storm impact regimes. The Grande Plage of Biarritz (GPB) is an urban embayed beach that is 1.2 km long, located on the southern Aquitanian coast of France (Figure 2). It has a high socio-economic importance for the city of Biarritz due to its tourist appeal, its historical heritage, and its location near the city center. In terms of characteristics, the GPB is an intermediatereflective beach with typically a steep foreshore slope of 8-9% and a gentle nearshore slope of 2-3% [START_REF] Morichon | Assessment of Flooding Hazards at an Engineered Beach during Extreme Events: Biarritz, SW France[END_REF]. It is a mesotidal beach with 4.5 m spring tidal range around a mean water level of 2.64 m. This narrow beach is backed by a seawall with an alongshore elevation varying between 7 and 8 m. This seawall serves as defense infrastructure for back beach buildings.

The beach is predominantly exposed to waves coming from the WNW direction. The offshore wave climate is moderately to highly energetic. The annual average significant wave height and peak period are, respectively, H s = 1.5 m and T p = 10 s [START_REF] Morichon | Assessment of Flooding Hazards at an Engineered Beach during Extreme Events: Biarritz, SW France[END_REF]. In this region, an event is qualified as a storm event when H s and T p are, respectively, greater than 3.5 m and 13.8 s. Such events correspond to 7.24% of the offshore wave climate [START_REF] Abadie | Wave Climate and Longshore Drift on the South Aquitaine Coast[END_REF] and are responsible for several overwash events each year.

This site has been equipped with a coastal video monitoring station since 2017. The station includes 4 cameras with different lenses to ensure the coverage of the entire beach with a sufficient spatial resolution. The cameras are operated by the open source software SIRENA [12]. For this site, one transect is monitored by the camera pointing to the beach and seawall location (transect Stack-Id01 in Figure 3). The timestack images correspond to pixel intensities recorded along this transect over 14 min with a sampling frequency of 1Hz. Among the 70,000 images of this database, only 8172 images were kept to be part of the ground truth dataset. Indeed, the timestacks generated in summer months were excluded as the human activities negatively affect the quality of the images. The images where the tide level was below 2.8 m were excluded as they corresponded to timestacks images without visible swash.

Timestack Image Preprocessing

The ground truth dataset was built by labeling the 8172 images. There are two methods to annotate the images: by hand or in a semi-automatic way. The annotation by hand is the most straightforward but also the most time-consuming method. The semi-automatic method consists of two steps. First, the position of the waterline is extracted automatically by segmenting the image using Otsu's thresholding method [16]. Then, the storm impact regime is identified by comparing the position of the waterline with the one of the defense infrastructure. This method is faster than the annotation by hand; however, it still requires an operator because it is not always robust and highly depends on the lightning conditions of the image. To employ this method, the position of the defense infrastructure in the image must be known. This is the case for the Grande Plage de Biarritz; therefore, semi-automatic annotation was performed.

After verification and correction by an operator, the result of the annotation was 7907/211/54 (Swash/Collision/Overwash). The classes are highly imbalanced, and this could have some effect on the classification accuracy of the CNN. Methods to deal with this problem are presented below. Before the training process, the images were resized to fit to the input dimensions of the CNNs tested in this study (224 × 224).

Zarautz 2.2.1. Site Characteristics

The beach of Zarautz is an embayed beach of 2.3 km long located on the Basque coast (northern Spain) in the SE Bay of Biscay, approximately 70 km southwest of GPB. The beach, facing north (345 degrees), can be divided into two parts (Figure 4): 30% of the beach in the eastern part presents a large and well-preserved dune system, with a maximum height of approximately 10 m above the minimum astronomic tide. The remaining 70% is an engineered urban beach, backed by a concrete seawall and the village of Zarautz.

In terms of characteristics, the beach of Zarautz is an intermediate-dissipative [START_REF] De Santiago | Video monitoring nearshore sandbar morphodynamics on a partially engineered embayed beach[END_REF] and mesotidal beach with a 4 m spring tidal range. It is composed of fine-medium sand with a mean slope of around 2%. The annual average significant wave height and peak period are, respectively, 1 m and 9 s. Like the GPB, the beach of Zarautz is also exposed to highly energetic waves and storms coming from the WNW and NW directions. The seawall backing 70% of the beach has an along-shore elevation varying between 6.5 m in the western part and 8 m in the center of the beach. This seawall serves as a defense infrastructure for the buildings near the beach, and overtopping events are common at high tide during winter storms.

A video monitoring station, like the one used on the GPB site was installed in 2010. The station has 4 cameras of 1.4 Megapixels. Two of the cameras are equipped with 12 mm lens and have a panoramic view, and another 2 equipped with 25 mm lens cover with more resolution the mean high and low tide coastline positions. For the Zarautz dataset, 4 transects are monitored by the camera covering the supra-tidal beach with higher resolution (Figure 4). The transects are perpendicular to the seawall and are named corresponding to the elevation of the seawall in the point of intersection (i.e., transect 65 corresponds to the part of the seawall with 6.5 m elevation).

Timestacks Images Preprocessing

Images from the site of Zarautz were annotated by hand. This method of annotation was preferred over semi-automatic annotation because: (i) the position of the seawall varied between timestack images, making the semi-automatic method more laborious, and (ii) the presence of strong winds and gust negatively impacted the quality of certain images, making the semi automatic method less robust. A simple web application was developed to facilitate the annotation for the operator and is accessible in a public GitHub repository (link in Data Availability Statement section). After classification by hand, the result of the annotation was 19,596/2776/162 (Swash/Collision/Overwash). Like the images of Biarritz, images of Zarautz present class imbalance, and they were resized to fit to the input dimensions of the CNNs before the training.

Convolutional Neural Networks

General Concept

CNNs are a type of neural networks widely used to perform tasks related to imagery analysis such as image segmentation, classification, or object detection. For classification problems, a CNN takes as input images with three channels (RGB), from which they output probabilities of belonging to specified categories, in our case storm impact regimes. Like a classical neural network, a CNN is a stacking of neurons that are organized in different layers. The structure of a CNN can be divided into two parts. The first part contains mostly convolutional and pooling layers and aims to learn specific features that help to classify the images correctly. The second part contains fully connected layers and the output layer. It uses the specific features extracted in the first part to output probabilities of belonging to specified categories.

In the feature extraction part, the convolutional layers detect features inside an image. They convolve their input with one or more filters, which results in one or more feature maps (one for each filter). The feature maps represent the activation of a specific filter at every spatial position of the input image. During the learning process, the network will learn filters that activate when they see specific visual features that help to correctly classify the training images. Usually, convolutional layers are stacked inside a CNN. The early layers detect simple features such as edges, whereas the deeper layers can detect more complex features.

Pooling layers are commonly found between convolutional layers. These layers also rely on convolutional operations and aim to reduce the dimensionality of the feature maps in order to increase the learning speed of the network and to control the overfitting of the CNN. If a CNN is overfitted, it would indicate that the network has learned exactly the characteristics of the training images and cannot generalize to new data. By stacking several convolutional and pooling layers inside a CNN, the complexity of the extracted information increases as we go deeper in the network with more feature maps with smaller dimensions.

The output of these specific layers serves as input to the second part of the network, which aims to classify the image into the correct category. In the classification part, neurons are organized in layers and are connected to the previous layers through weights (hence the name fully connected layers). To prevent overfitting, drop-out regularization can be applied on these layers. This method randomly ignores neurons during the training process, making the network learn more sparse and robust representation of the data. Finally, the output layer estimates the probabilities of belonging to the specified categories for the input image with a "softmax" activation function.

The CNNs are trained with backpropagation in the same manner as classical neural networks: the weights in the convolutional and fully connected layers are updated iteratively to minimize the errors between the prediction of the network and the ground truth. The ground truth dataset for such a network is made by annotating images. Details on the annotation of the timestacks can be found in the "Study Sites" section. Only the general ideas about CNN have been presented above; for a detailed description on CNN and their training, the reader is referred to the work of Bengio et al. [START_REF] Bengio | Deep Learning[END_REF].

There are many CNN architectures, each with different complexity and characteristics. In order to keep the computation time reasonable, it was decided to limit the comparison of performances between four architectures of increasing depth and complexity:

•

A custom architecture inspired by the work of LeCun et al. [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] adapted for bigger images. The architecture is presented in the appendix (Table A1). • AlexNet [START_REF] Krizhevsky | Imagenet Classification with Deep Convolutional Neural Networks[END_REF], which won the ImageNet challenge in 2012. Its architecture contains more convolutional layers and dense layers (Table A2). The number of filters is also larger than that of the custom architecture. • VGG16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF], which is a very deep CNN that uses 13 convolutional layers and three dense layers (Table A3). • Inception v3, an improved version of the GoogleNet from Szegedy et al. [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF] which won the ILSVRC in 2014. It relies on inception modules, which perform convolutions with filters of multiple size and concatenate their results (Table A4). In addition, the convolution operation with filters of large size inside an inception module are made by using 1 × n filters to reduce computational cost. This results in deeper networks with significantly fewer parameters to learn.

Training the CNN

Data Processing

The datasets of both sites were divided into training, validation and testing sets containing, respectively, 65%, 15%, and 20% of the data (common proportions in the literature). Stratified random sampling was used to ensure that each part contains the same class proportions. The training set is used to fit the CNN. The validation set is used to stop the training for the CNN (early stopping). At last, the test part is used to evaluate the performance of the neural network on unseen data (not used in the training step).

During the training, each training image is seen multiple times by the CNN. This can be a problem as the network can learn exactly the characteristics of the training images and might not generalize to new data. To avoid this problem, called overfitting, data augmentation is employed during the training of the CNN. This method consists in making small changes to images in the training set before feeding them to the CNN. By generating modified images, this method artificially increases the number of images in the minority classes and makes the models more robust to overfitting.

The following changes were made: The datasets from both sites suffer from the class imbalance problem. Indeed, the distributions of storm impact regimes are highly imbalanced. For the Biarritz site, 96.8% of the images display swash regimes, 2.6% collision regimes, and 1% of overwash regimes. For Zarautz, 87% of the images are swash regimes, 12% collision regimes, and 1% overwash regimes. This class imbalance problem was expected as we are studying rare events.

It has been proven that class imbalance can negatively affect the performances of machine learning models in classification tasks [START_REF] Japkowicz | The class imbalance problem: A systematic study[END_REF]. Methods to deal with this problem are well known [START_REF] Japkowicz | The class imbalance problem: A systematic study[END_REF][START_REF] Buda | A Systematic Study of the Class Imbalance Problem in Convolutional Neural Networks[END_REF][START_REF] Johnson | Survey on deep learning with class imbalance[END_REF] and can be divided into two categories: data-level methods and classifier-level methods.

The data-level methods aim to modify the class distribution in order to reduce the imbalance ratio between classes. The most popular methods in this category are oversampling and undersampling. Oversampling consists in replicating random samples from minority classes until the imbalance problem is removed. In contrast, undersampling consists in removing random samples from the majority class until the balance between classes is reached.

The classifier-level methods aim to modify the training or the output of the machine learning algorithm. They include cost-sensitive learning, which is a method that gives more weights during learning to examples belonging to minority classes, and the thresholding method, which adjusts the output probabilities by taking into account the prior class probabilities [START_REF] Buda | A Systematic Study of the Class Imbalance Problem in Convolutional Neural Networks[END_REF].

Transfer Learning

For complex and deep CNN, it is common to use transfer learning to speed up the learning process and to improve performances. Transfer learning methods consist in using knowledge gained on a specific task to solve a different task. There are different methods of transfer learning for CNN; for an exhaustive listing, readers are referred to the work of Pan and Yang [28]. The method used in this article is "pre-training". It consists of using the weights of a CNN trained on a first task as initialization weights for a second CNN that will perform on a second task. The efficiency of pre-training was tested by using the pre-trained weights of VGG16 and Inceptionv3 on the ImageNet dataset, which is one of the largest labeled image dataset [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF]. Then, transfer learning was performed between sites to see if the knowledge gained on one site is beneficial for the learning on the second site.

Application to the Datasets

The workflow for this study is presented in Figure 5. For each site, the four CNNs with different architectures were fitted without and with the two methods related to the class imbalance problem: oversampling and cost-sensitive learning (class weights). Transfer learning was used on the more complex architectures (VGG16, Inceptionv3) and only for the best performing method to cope with class imbalance. Data augmentation was used during the training of all the CNNs.

The networks were trained on a laptop equipped with a GPU (Quadro RTX 4000) using Keras (tensorflow GPU 1.12.0/Keras 2.3.1/Python 3.6. 1), an open-source python library designed for building and training neural networks. The scripts used in this article are available on a public GitHub repository (link in Data Availability Statement section). The optimizer used is Mini-Batch gradient descent algorithm with batch size of 32 and a learning rate of 0.001 that decays by a factor of 2 every 10 epochs. The training is stopped at 100 epochs or earlier when the value of the validation loss does not decrease over 10 epochs (early stopping).

CNN Accuracy Assessment

To compare the performance of the different networks, the F 1 -score is computed with the following formula: The precision, recall, and F 1 -score are computed for each storm impact regime and are averaged in order to have one global metric for each CNN. The F 1 -score varies between 0 and 1, with 1 representing the best value. Unlike the global accuracy (number of correct predictions divided by the total number of predictions), the F 1 metric is not biased when data present a class imbalance.

Results

The results are organized into four subsections. Firstly, the performances of the different combinations of CNN architectures, methods to cope with class imbalance, and transfer learning are compared. Secondly, the prediction errors of the best CNN for each site are investigated. Thirdly, we present results related to transferability between sites. Finally, a sensitivity analysis is presented for the site of Zarautz.

CNN Performances

Architectures

Table 1 regroups the training time, number of epochs, and also the performance metrics (accuracy, recall, F 1 -score) for different CNN architectures, methods to cope with class imbalance, and with or without pre-training with ImageNet dataset for both sites. For both sites and for every methods used to cope with class imbalance (or not), CNNs with deeper and more complex architectures yielded better results (higher values of precision, recall, and F 1 -score). Indeed, these kind of architectures tend to learn more complex features that lead to better performance in harder tasks. The downside of these networks is the training time, which is significantly higher than simpler and shallower models. Without coping with class imbalance problem, CNNs tend to predict all the images as the majority class, resulting in poor classification results. Between the two methods tested, oversampling seems to perform better, with F 1 -scores on average 30% better than the ones obtained with cost-sensitive learning (class weights). The superior performance of oversampling method on this dataset might be due to the fact that the CNNs see more images during training in the oversampling method than in the cost sensitive learning method, resulting in better classification accuracy. Finally, models using pre-trained weights (transfer learning) train faster (fewer epochs) and yield better classification results than models trained from scratch. Indeed, the F1-scores obtained with the pre-trained models are, respectively, 6 to 8% higher for the VGG16 and Inceptionv3 models. Even though the images from the ImageNet dataset have different characteristics than the timestack images that are being classified, the pre-trained weights might contain knowledge about general features that are helpful to better classify the timestacks.

Best Models

For GPB, the best model is the pre-trained VGG16 with an F1-score of 0.866. The pretrained Inception model trains faster but shows a lower F1-score (0.780). For the Zarautz site, the best model is also the pre-trained VGG16 with an F1-score of 0.858, but this time the performance of the pre-trained Inception v3 model was very close with an F1-score of 0.849.

Investigating the Errors

The confusion matrices on the test sets are presented in Table 2. In general, the minority classes tended to have higher error rate. This is expected as minority classes contain fewer examples than majority classes. The prediction errors made by the CNN were manually inspected to gain an understanding of common error types. Prediction errors made on the GPB test set are presented in Figure 6 and in the appendix (Table A5). Among the 16 errors made on the test set, five errors came from human misclassification, five errors may have been caused by the presence of specific features such as vertical lines usually associated with collision regimes, two errors were made on images that were displayed in between category of storm impact regimes. The remaining errors correspond to images that were hard to classify, where lighting and meteorological conditions were poor. The misclassification errors were corrected for the test and validation sets, and the best network was trained once again, resulting in slightly lower results but this time without human misclassification error (Table 3). The errors made on the Zarautz dataset were also analyzed (Table A6). A large number of errors were made on images that were in between the categories of storm impact regimes: either the images displayed a swash regime, which was very close to the collision regime, or the images displayed a regime impact, with one small overtopping of the wall. Some misclassification errors were made. The rest of the errors may have come from lighting conditions (large horizontal band, lighter in the images). The misclassification errors were corrected for the test and validation, and the best network was trained once again, resulting in a slightly better result for this site (Table 3).

Transferability between Sites

The interest in transfer learning for CNN training has been highlighted in Section 4.1.3. The models using pre-trained weights from ImageNet trained faster (fewer epochs) and yielded better classification results. In this section, we investigate if the knowledge acquired on one site (CNN weights) could be transferred to another site with pre-training. Pretraining is the most common way to transfer knowledge between tasks. It consists in using the weights of a neural network trained on one site as initialization weights for the training on the second site. The weights of the best CNN for each site are used for the pre-training of a CNN on the other site. The performances of these CNN are presented in Table 4. The weights of the best model on Zarautz data (VGG16 transfer) were used as initialization weights for the training on Biarritz data. This resulted in classification results better than the learning from scratch with a higher precision and F 1 -score (Table 4). However, the values of precision, recall, and F 1 -score obtained with pre-training on Zarautz data remained slightly lower than the ones obtained with pre-training on ImageNet data.

Pre-training method was also applied on Zarautz data, where the weights of the best model on the Biarritz site were used as initialization weights. The classification results were better than learning from scratch and learning with pre-trained weights from ImageNet with higher F 1 -score (Table 4).

Sensitivity Analysis

A sensitivity analysis was performed on the dataset of Zarautz to highlight the effect of the size of the training images dataset on the classification accuracy. The dataset of Zarautz was divided into three smaller datasets. Each of these datasets was divided into the training/validation/test sets with the proportions described in Section 3.2.1. Finally, a CNN model was trained on each smaller datasets (VGG16 with transfer learning ImageNet and oversampling).

The averaged F1-metric for these three models was 0.805. This value is slightly lower than the one obtained with the full dataset, which is 0.858 (Table 1). These results confirm what was already known in the literature: CNN performances tend to increase with the training set size [START_REF] Uchida | A further step to perfect accuracy by training CNN with larger data[END_REF].

Discussion

Even though we showed the strong potential of CNN to automatically generate storm impact regime database, the proposed methodology can be improved in several ways. More attention could be paid to the choice of CNN architectures and hyperparameters. Other CNN architectures need to be tested, especially recent architectures such as ResNet, MobileNet, or Xception. They could perform better than the architectures presented in this work. For instance, the ResNet architecture contains skip connections between layers, which allows the training of much deeper and performant networks. Hyperparameters are parameters whose values are specified by the user before the training process begins; they affect the structure of a CNN and how well it trains. They have a non negligible impact on the final results. Several optimization algorithms such as Bayesian optimization could be employed to select the optimal hyperparameters [START_REF] Snoek | Practical Bayesian Optimization of Machine Learning Algorithms[END_REF], which has not been done in this study.

In addition to hyperparameter tuning, other methods of data augmentation could be used to improve the performances of the CNN. The analysis of prediction errors can help in the choice of other data augmentation methods. In our case, many errors were related to lighting conditions; it would be wise to test various data augmentation methods affecting the lighting or brightness of the images. This could make the CNN more robust to lighting conditions and therefore improve its performances.

It is worth noting that the performances of a CNN model implemented at a given site are expected to increase with time as more timestacks are collected by the video monitoring system. With more training images, the minority classes will contain more images, and this will lead to less classification errors for these classes. Moreover, if enough timestacks are collected, intermediate storm impact regime classes could be created. These classes could reduce the errors on the images displaying impact regimes not corresponding to the three regimes presented in this work.

One very interesting feature of CNNs models is their transferability. We showed that using the knowledge acquired from another site can lead to improved classification results when using pre-training (especially for Zarautz site). The weights of the best CNNs for both sites are available in a GitHub repository (link in Data Availability Statement section) and could be used as initialization weights for a CNN applied to a new site. The only requirement is to annotate timestacks from the new site, which will serve as training data.

Despite the promising performance, this methodology has some limitations, mainly related to the image annotation, an obligatory step for CNN training. The first limitation of this method is the lack of knowledge about its sensitivity. We showed for the site of Zarautz that CNNs yield lower performances when trained on a smaller training set. However, we do not know the minimum number of timestacks to annotate for a new site in order to have satisfactory accuracy. A sensitivity analysis should be performed to find this minimum threshold and to make some recommendations on the use of this method in the case of new sites with a small number of timestacks.

The second limitation of this method is the annotation process itself, which is tedious and time-consuming. An alternative solution could be to use the domain-adaptation approach presented in the work of Ganin and Lempitsky [START_REF] Ganin | Unsupervised domain adaptation by backpropagation[END_REF]. They propose a specific CNN architecture that can be trained simultaneously on a large number of labeled data from a source domain (one site) and unlabeled data from a target domain (new site). At the end of the training, the CNN is able to classify correctly images from both sites even though only images from one site have been labeled.

Finally, the performances of the proposed method must be compared objectively with human-level performance and other methodologies. Assessing the human-level performance on this task is essential and would give precious insights into how to improve further the CNN performances [START_REF] Dodge | A study and comparison of human and deep learning recognition performance under visual distortions[END_REF]. For example, a CNN performance lower than the human-level performance could indicate the presence of a bias, which can be avoided by using deeper models or by training more slowly and for longer. It would be of great interest to compare this methodology based on CNNs with methodologies based only on traditional imagery analysis. As stated in the introduction, a possible methodology could be to first extract the waterline position using Otsu's segmentation [16,26] or using the radon transform [27] and then compare its position with the position of defense infrastructure to define the storm impact. Another methodology could be based on the analysis of pixel intensity such as the works of Simarro et al. [START_REF] Simarro | On the use of variance images for runup and shoreline detection[END_REF] and Andriolo et al. [START_REF] Andriolo | Operational use of surfcam online streaming images for coastal morphodynamic studies[END_REF]. The methodologies based on simple image processing algorithms could have some advantages over CNNs. Indeed, they would not require the building and training of a CNN structure, which is time-consuming, and the whole decisional process is known to be contrary to CNN, which can be considered as a "black box". However, these methodologies would need to be adapted for each site by indicating the position of the defense infrastructure, which is not needed with CNN. In addition, the simple image processing algorithms could be more affected by the erratic brightness of the timestacks than CNNs, which are trained with data augmentation.

This work is a first step in the analysis of storm impact with video monitoring. Numerous extensions can be envisaged, particularly on the type of information extracted and the type of image analyzed. Indeed, the CNN could be used to count the number of collision or overwash events in one timestack. This technique could be also extended to analyze other types of images produced by video monitoring systems such as oblique and/or rectified images. Finally, it can be employed to analyze images from already existing cameras such as surfcam [START_REF] Mole | Capitalizing on the surfcam phenomenon: A pilot study in regional-scale shoreline and inshore wave monitoring utilizing existing camera infrastructure[END_REF][START_REF] Andriolo | Nearshore wave transformation domains from video imagery[END_REF]. This could constitute a low-cost monitoring method with a large spatial coverage for the qualitative study of storm impact. Many questions arose with this work, especially about the minimum number of images to annotate to have satisfactory accuracy or the lack of comparison with the current method or human level performance. More questions will arise during the operational implementation and use of the CNNs concerning the verification of predictions, the prediction error handling, or how often we need to re-train the neural networks with the newly classified images.

Conclusions

In this paper, we presented an innovative methodology based on convolutional neural networks and coastal imagery that could be used to collect storm impact data routinely. We described the methodologies associated with CNNs, including the annotation of the dataset, the training of the networks, or transfer learning. We also introduced the problem of class imbalance, which is due to the extreme nature of the storm impact regimes, and we proposed and compared different solutions such as oversampling or cost-sensitive learning.

The proposed methodology was tested on two sites: Biarritz and Zarautz. We showed that convolutional neural networks are well adapted for the classification of timestacks into storm impact regimes. Overall, we found that more complex and deeper architectures yielded better results. Best performances were achieved with the VGG16 architecture for both sites with F-scores of 0.866 for the site of Biarritz and 0.858 for the site Zarautz. For the class imbalance problem, the method of oversampling showed better classification accuracy than the cost-sensitive learning method, with F-scores on average 30% higher. Finally, we showed that the method can be easily applied to a new site with optimal efficiency using transfer learning. Indeed, training a CNN using pre-trained weights (ImageNet or weights of another site) resulted in better accuracy than training a CNN from scratch (F-scores on average 6 to 8% higher).

With convolutional neural networks, we can take full advantage of the large number of data produced by video monitoring systems. We showed that they are able to transform images into usable qualitative data about storm impact. Even if the data are not continuous (only day time and winter months), this method could be, without a doubt, a real asset in the future for coastal researchers and stakeholders by routinely collecting storm impact data, which are rare at present. These data are essential in the disaster risk reduction chain, and they have many uses. They can serve as validation data for impact models or early warning systems based on numerical modeling. They can also be used to train early warning system based on Bayesian networks [START_REF] Poelhekke | Predicting Coastal Hazards for Sandy Coasts with a Bayesian Network Coast[END_REF][START_REF] Plomaritis | Use of a Bayesian Network for Coastal Hazards, Impact and Disaster Risk Reduction Assessment at a Coastal Barrier (Ria Formosa, Portugal)[END_REF]. Finally, statistical analysis can be performed to find relationships between observed storm impact regimes and local conditions such as wave characteristics, tide, or meteorological conditions. funding acquisition, D.M. and B.L. All authors have read and agreed to the published version of the manuscript.

•

Total params: 14,978,883 • Trainable params: 14,978,883 • Non-trainable params: 0 Table A6. Errors explanation for Zarautz data. "Misclass." stands for misclassification during annotation, "Splash" corresponds to an intermediate storm regime between impact and overwash, "Vertical" corresponds to the presences of vertical features of runup. Finally, "SI" corresponds to an intermediate storm regime between swash and impact that was very close to the sea wall but did not impact.

Conclusion

This chapter showed the ability of deep learning methods, in particular Convolutional neural networks (CNN) to create automatically a storm impact database with images from video monitoring stations. We demonstrated that deeper CNN architectures associated with oversampling method yielded best classification results for two study sites which are Biarritz and Zarautz. This methodology is easily transferable on new sites by using the pre-trained weights of other sites. Even though the database constituted by this methodology is only qualitative, it can already be used to train or validate statistical models aiming to predict storm impact. In the future, CNN could be essential assets as the monitoring networks are expected to be more dense with new technological advances and lower cost of the components [START_REF] Xu | Internet of things in marine environment monitoring: A review[END_REF][START_REF] Marcelli | Toward the widespread application of low-cost technologies in coastal ocean observing (internet of things for the ocean)[END_REF].

Many extensions of this work can be envisaged. As of now, this methodology only extracts maximum storm impact regime for a timestack image resuming 15 minutes of video. A lot of information about coastal flooding are therefore not taken into account such as the number of collision or overwash events, or the time when they happen. One possible extension of this work could be to detect and count each collision/overwash in each timestack, leading to better risk characterization.

Scientific output:

This work resulted in a communication during an international conference:

• "Developing a Bayesian network to predict coastal flooding on the Grande Bayesian networks to model storm impact using data from both monitoring networks and statistical learning methods

Abstract

Bayesian networks (BNs) are probabilistic graphical models that are increasingly used to translate hydraulic boundary conditions during storm events into onshore hazard. However, comprehensive databases representative of the extreme and episodic nature of storms are needed to train the BNs. Such databases do not exist for many sites and many BNs are trained on data generated by process-based model. To our knowledge, BNs have not been trained exclusively on observational data in coastal engineering domain.

This study aims to explore the performance of a BN exclusively based on observational data in coastal flooding prediction. To this end, we use as a test case the Grande Plage of Biarritz (South west of France). The BN is trained using data from several monitoring networks located near the study site. Because observational data about storm impact regime and atmospheric surge are limited, a second aim of this work is to propose a methodology based on statistical learning methods to extend the data about these variables. This methodology aims to select the statistical learning method with the best generalizing ability with a cross validation. Two BNs are trained, one exclusively on the observational data and one with both observational and predicted data. To compare the two networks, their performances are evaluated on the same events.

We demonstrated that it was possible to predict coastal flooding risk in a qualitative manner with a BN based only on observational data with a F 1 -score of 0.628. However, the predictive skill of this network is questionable for the most intense storm impact regimes which are impact and overwash regime. Storm impact and atmospheric surge data were both extended by random forest method which is the method that showed the best generalizing ability in the two cross validation. This extension of the database led to a better BN in terms of predictive skill, with precision, recall and F 1 -score in average 20% higher than the BN

Introduction

The frequency and intensity of coastal flooding is expected to increase in the future due to sea level rise and the continuous development of coastal areas (Vousdoukas et al., 2018;Taherkhani et al., 2020). The development of early warning systems (EWS) is therefore mandatory to mitigate the risk related to coastal flooding events. A crucial step in an operational EWS is the translation of the hydraulic boundary conditions into onshore hazard. This translation is performed most of the time by phase-resolving wave models. This type of wave models is employed to describe the free surface elevation at the scale of single waves. They have the ability to represent complex processes such as wave setup, re-circulation, and infra-gravity waves as a result of the free surface evolution [START_REF] Roeber | Improving wave run-up forecasts-benefits from phase-resolving models[END_REF]. However, these models require a fine resolution to perform well and are therefore computationally expensive. This complicates their use in EWS which requires fast predictions. Recently, more attention have been given to bayesian networks (BNs) for the development of the storm impact model.

BNs are a class of probabilistic graphical models which allow for an intuitive representation of a set of random variables and their conditional dependencies [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF][START_REF] Scutari | Bayesian networks: With examples in r[END_REF]. They can learn causal effect from observational data, represent complex systems with intuitive graphical structure and include uncertainty. For these reasons, applications of BNs have multiplied recently in various domains spanning from ecology [START_REF] Mccann | Bayesian belief networks: applications in ecology and natural resource management[END_REF][START_REF] Landuyt | A review of bayesian belief networks in ecosystem service modelling[END_REF] to coastal engineering (Poelhekke et al., 2016;[START_REF] Jäger | A Bayesian network approach for coastal risk analysis and decision making[END_REF]Beuzen et al., 2018).

In coastal engineering, BNs are appreciated methods due to their low computational cost and their ability to represent complex systems by integrating different sources of data. Their intuitive representation is also a non negligible advantage compared to other modeling approaches (e.g. process-based models) as it facilitates the communication between scientists and stakeholders in management applications [START_REF] Henriksen | Public participation modelling using bayesian networks in management of groundwater contamination[END_REF]. Most of the applications of BNs in coastal engineering domain concern the translation of forcing variables (e.g. wave, weather and tide conditions) into impact and damages on the shore during storm events. Indeed, they have been employed to predict coastal cliff erosion (Hapke and Plant, 2010), shoreline retreat (Beuzen et al., 2018), dune retreat and erosion (Palmsten et al., 2014;den Heijer et al., 2012) and barrier island response (Plant and Stockdon, 2012;Wilson et al., 2015;Poelhekke et al., 2016) resulting from coastal storms.

A challenge in the development of BNs for storm impact modeling is related to the training database. The training database must be comprehensive and cover a long time period in order to be representative of the extreme and episodic nature of storms. However, such databases do not exist for many coastal sites, therefore, previous applications of BNs rely on process-based modeling to either simulate the forcing parameters or to translate the forcing conditions into impact on the shore using a numerical model to simulate nearshore waves transformation and their impact at the coast. For instance, Poelhekke et al. (2016) and Plomaritis et al. (2018) create synthetic storm events by using copula statistical method fitted on the limited observations of offshore hydrodynamic parameters. The impact of these synthetic events are then simulated using the nearshore wave propagation model Xbeach [START_REF] Roelvink | Xbeach model description and manual[END_REF]. The BNs are then trained to link the characteristics of the synthetic storm events to their simulated impact. Lately, with the increasing amount of data collected by diverse monitoring networks, observational data are being incorporated in the training data of the BNs. In Beuzen et al. (2018), the BN is trained on 10 years of storm wave data simulated by process-based models (spectral wind wave model) and also on shoreline retreat data derived from a coastal imaging station.

To our knowledge, BNs have not been trained exclusively on observational data in coastal engineering domain. This study aims to explore the performance of a BN exclusively based on observational data in coastal flooding prediction. To this end, we use as a test case the Grande Plage of Biarritz (South west of France). Every winter, this embayed beach faces numerous storm events which often result in coastal flooding. To prevent and mitigate the risk of flooding, a BN is developed using data from several monitoring networks located near the study site. Forcing parameters are collected by a directional wave buoy, a tide gauge and a weather station. Data about storm impact on the beach are determined from the video monitoring station installed on the site, following the methodology presented in [START_REF] Callens | Automatic creation of storm impact database based on video monitoring and convolutional neural networks[END_REF].

Due to the recent installation of the video monitoring system (in 2017), the storm impact data is limited. It is why we propose a methodology based on machine learning methods to extend the database. The potential gain of performance with the extension of the database is investigated by building two different BNs: one only using observational data and the other based on both observational data and data predicted by the ML algorithms. Their performances are compared by testing their accuracy on the same events. Section 2 will introduce the study area, section 3 will introduce the theory behind bayesian networks and the data used to train them. Results will be presented and discussed in Section 4. Finally, Section 5 will cover the conclusion.

Study site

The Grande Plage of Biarritz (GPB) is an urban embayed beach of 1,2 km long located on the southern Aquitanian coast of France (Figure 1). It has a high socio-economic importance for the city of Biarritz due to its tourist appeal, its historical heritage and its location near the city center. In terms of characteristics, the GPB is an intermediate-reflective beach with typically a steep foreshore slope of 8 -9% and a gentle nearshore slope of 2 -3%. It is a mesotidal beach with 4.5 m spring tidal range around a mean water level of 2.64 m. This narrow beach is backed by a seawall with an alongshore elevation varying between 7 and 8m. This seawall serves as defense infrastructure for back beach buildings. The beach is predominantly exposed to wave coming from the WNW direction. The annual average significant wave height and peak period are respectively H s = 1.5m and T p = 10s. Storms in the region are considered as those in Hs and Tp are respectively greater than 3.5 m and 13.8 s. Such events correspond to 7.24% of the offshore wave climate (Abadie et al., 2006) and are responsible for several coastal flooding events each year.

This site is equipped with a coastal video monitoring station since 2017. The station includes 4 cameras with different lens to ensure the coverage of the entire beach with a sufficient spatial resolution. The cameras are operated by the open source software SIRENA (Nieto et al., 2010).

Methodology

Bayesian networks

In this section, we give a general overview of the theory behind BNs. The reader is referred to [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]; [START_REF] Scutari | Bayesian networks: With examples in r[END_REF] for in-depth details and examples of applications of BNs .

A bayesian network represents a set of random variables X = {X 1 , ..., X n } and their dependencies via a directed acyclic graph [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. In this directed acyclic graph (DAG), each variable is represented as a node. The nodes can be connected together by arcs that represent potential dependence between the variables. This arc is directed depending the direction of the influence from parent to child node. In a BN, the arcs must not form a cycle, hence the name of directed acyclic graph.

A fundamental property of the BN is that the global distribution of the set of random variables can be economically factorized through the chain rule [START_REF] Jäger | A Bayesian network approach for coastal risk analysis and decision making[END_REF]:

P(X 1 , ..., X n ) = n i=1 P(X i |pa(X i )), (1)
with pa(X i ) denoting the parent nodes of X i and P(X i |pa(X i )) specified by the conditional probability tables (CPTs), or the probability tables (PTs) when the nodes do not have parents. In discrete bayesian networks, the variables are discretized into states, in this case, the CPTs and PTs quantify the probability of a node being in its particular states.

Given an evidence, the BN can update the states of all the nodes connected to this evidence. This update relies on Bayes' rule:

p(R i |O j ) = p(O j |R i )p(R i ) P(O j ) (2)
where p(R i |O j ) is the updated or conditional probability of the response R i given O j a set of observations/evidences, p(R i ) is the probability of response R i happening, p(O j ) represents the probability of the observations O j .

The structure, the CPTs and PTs of a BN can be learned from a dataset. In this work, the structure of the BNs and their parameters are learned with the academic version of GeNIe Modeler software from BayesFusion (http://www. bayesfusion.com/). In this software, the structure of the BNs is learned using Bayesian Search which is a score-based algorithm [START_REF] Bayesfusion | Genie modeler. User Manual[END_REF]. This learning algorithm is based on hill-climbing and uses the out-of-sample classification accuracy (computed with 5-fold cross validation) as the scoring function to find the optimal graph. For parameters estimation (learning CPTs and PTs from the data), GeNIe Modeler uses the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF].

Data sources

Table 1 summarizes the characteristics of the data available for the study site. The majority of the variables has an hourly frequency. However this is not the case for the buoy measurements and the storm impact regime variable. Each data source is described more precisely below.

Wave data

We are interested only in the three wave integrated parameters that play a major role in coastal flooding: the significant wave height (H s ), the peak period (T p ) and the peak wave direction (θ p ).

Measured data

Direct measurements of these parameters are obtained from the National Center for Archiving Swell Measurements (L'her et al., 1999). They were made by a directional wave rider buoy (DWR MKIII) operated by the Centre for Studies and Expertise on Risks, Environment, Mobility, and Urban and Country Planning (CEREMA). The buoy is located a few miles off the Basque Coast (Figure 1) at 50 meters water depth. Since its deployment in 2009, this buoy have been recording the parameters of interest every 30 minutes.

Reanalysis data improved with error prediction method

The wave parameters are also simulated at the buoy coordinates using two hindcast sea-states database computed with the spectral wave model MFWAM developed by MeteoFrance (Lefèvre and Aouf, 2012), each covering two distinct period between 1993 and 2021:

• "IBI MULTIYEAR WAV 005 006". This reanalysis covers the period 1993-2019 with a hourly time-step. It is forced by the ERA 5 reanalysis wind data from ECMWF.

• "IBI ANALYSIS FORECAST WAV 005 005" covers the period 2020-2021.

It is forced with the ECMWF hourly wind data (forecast).

The spectral wave models are known to underestimate wave characteristics during energetic conditions (Rakha et al., 2007;Moeini et al., 2012;Arnoux et al., 2018). When measured data are available, it is common to perform data assimilation to improve the values of the wave parameters. The error prediction method presented in Callens et al. (2020) is therefore applied to both models in order to improve the predicted wave parameters. The performances before and after the data assimilation on both reanalysis are shown in appendix (Tables A1 andA2).

Tide data

The water level data come from the tide gauge located a few kilometers south of the study site (Figure 1) and were provided by the french Naval Hydrographic and Oceanographic Service (SHOM). The data range from 2011 to nowadays with a hourly time step. In addition to provide the data, the SHOM also documented each dysfunction of the tide gauge. This allowed for the proper removal of aberrant data. From tide data measurements, the astronomical tide was estimated with harmonic analysis performed by the R package TideHarmonics. The atmospheric surge was calculated by subtracting the astronomical tide to the measured water level.

Meteorological data

Meteorological data, including average wind speed above 10 meters, wind direction and atmospheric pressure were furnished by the French national meteorological service MeteoFrance. The data were collected hourly by the meteorological station of the Biarritz airport, located only a few kilometers from the study site (Figure 1). It covers a period ranging from 1993-01-01 to 2020-10-31.

Storm impact data

The storm impact regimes for the study site of Biarritz have been extracted from timestack images created by the video monitoring system of the Grande Plage of Biarritz (Figure 1). They are extracted with convolutional neural networks (CNN) following the methodology presented in [START_REF] Callens | Automatic creation of storm impact database based on video monitoring and convolutional neural networks[END_REF]. In this method, each timestack image is classified into 3 storm impact regimes representing increasing categories of coastal flooding risk:

• Swash regime: all the waves are confined to the beach • Impact regime: at least one wave collides with the bottom of the seawall • Overwash regime: at least one wave completely overtops the seawall The storm impact data ranges from 2017-03-23 to 2021-03-28. It contains 9550 swash regimes, 220 impact regimes and 54 overwash regimes.

Statistical models for database extension

Data about storm impact and atmospheric surge cover a shorter period compared to the other variables. In order to extend the data of these two variables, statistical learning methods (SLMs) are trained on the available data and employed to predict the occurrences where observations are missing.

For each variable, the performances of several SLMs are compared through a 5-fold cross validation on the available data. The aim of this comparison is to find the method with the best generalization performance. The hyperparameters of these methods are calibrated with bayesian optimization with the aim to maximize the generalization performance. For each variable, the model (and associated hyperparameters) showing the best generalization performance is used to predict the data when observations are not available. The cross validation and hyperparameter search is performed with the R package Tidymodels. Once the best predictive model is found, it is trained on all the available data and can be employed to predict data on the period of interest.

Concerning the atmospheric surge model, the explanatory variables are: the astronomical tide, the meteorological conditions (atmospheric pressure, wind speed and direction) and wave characteristics improved by the data assimilation method (H s , T p , θ p ). Three statistical learning methods are compared: gradient boosting trees, random forest and shallow neural networks.

For storm impact model, the explanatory variables are the astronomical tide, the wave characteristics improved by the data assimilation method (H s , T p , θ p ) and the atmospheric surge. Once again, three SLMs are compared through cross validation: gradient boosting trees, random forest and multinomial models.

Training the BNs

In order to assess the gain of performance with the extension of the database with statistical learning methods, we have to build two bayesian networks: one only based on observation data (Figure 2a) and the other based on observation and predicted data (Figure 2b).

The observational database is built by gathering the wave characteristics at the buoy, the observed atmospheric surge and astronomical tide extracted from the tide gauge, the weather conditions and the storm impact extracted from video monitoring network (Figure 2a). In total, we have 6358 observations with 6166 swash events, 151 impact and 41 overwash events. The test set is built with random stratified sampling and represents 20% of the observations (1234 swash, 31 impact and 9 overwash). The training set contains 4932 swash events, 120 impact and 32 overwash events.

The database extended with statistical learning methods is built by gathering the wave characteristics from the improved wave predictions, the astronomical tide extracted from the tide gauge, the weather conditions, the atmospheric surge and storm impact data (Figure 2b). For the the atmospheric surge and storm impact data, predictions of the statistical learning methods are included when observations are not available. The events of the test set are removed from the extended database to ensure an objective performance evaluation. The training set contains 237479 swash events, 860 impact and 292 overwash events.

Before training the bayesian networks, all the variables are discretized. Because the methods available to automatically discretize variables often find meaningless or questionable thresholds [START_REF] Chen | Predicting typhoon-induced storm surge tide with a two-dimensional hydrodynamic model and artificial neural network model[END_REF], the discretization of the variables in this study is based on expert judgment and visible thresholds observed in the data. The discretization thresholds for each parameter can be found in appendix (Table A5). In addition, partial undersampling is performed on both training set to minimize the imbalance ratio between class. Partial undersampling consists in randomly sampling a percentage of the observations belonging to the majority class.

It is an effective method to deal with the negative effect of class imbalance on the training of statistical learning methods (Japkowicz and Stephen, 2002). For the observational training set, only 20% of the swash events are kept resulting in 86.7% swash events, 10.5% of impact events and 2.8% overwash events. For the training set containing both observational and predicted data, only 2.5% of the swash events are kept to match the class distribution of the first training set (83.8% swash, 12.1% of impact and 4.1% overwash).

Evaluating the BNs

The descriptive and predictive skill of both BNs are evaluated. The descriptive skill corresponds to the predictive performance of the BN on the training set while the predictive skill corresponds to the predictive performance on unseen data (test set). The performances are evaluated with F 1 -score, precision and recall: These metrics are computed for each storm impact regime and are averaged in order to have one global metric for each BN. The F 1 -score varies between 0 and 1, with 1 representing the best value. Unlike the global accuracy (number of correct predictions divided by the total number of predictions), the F 1 metric is not biased when data present a class imbalance.

Results

Statistical models to extend the database

Atmospheric surge model The performances of random forest, gradient boosting trees and neural networks are presented in Table 2. The optimal hyperparameters found for these algorithms are shown in appendix (Table A3).Random forest shows the best performance for atmospheric surge modeling with a mean RMSE of 0.0565. This method is followed by gradient boosting trees and neural networks. Random forest with the optimal hyperparameters is trained on the totality of the atmospheric surge observations (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) and then used to predict the period 1993-2011 that does not have observations for atmospheric surge.

Storm impact model

The performances of the three SLMs for storm impact prediction are presented in Table 3. The optimal hyperparameters found for these algorithms are shown in appendix (Table A4) Kappa metric was used to compare the performances of the SLMs. This metric measures the degree of agreement between the true values and the predicted values by a classifier and has the advantage to not being affected by class imbalance. A value of 1 represents a perfect agreement between predictions and true values whereas 0 represents chance agreement. Here again, random forest showed the best generalization performance with a kappa metric of 0.68. Random forest with its optimal hyperparameters is trained on all the storm impact data available (2017-2021). To have an idea about the descriptive performance of this model, it has been used to predict on the training set. The confusion matrix between the predictions of the model and the real observations of storm impact is presented in Table 4. The predictions of this model are not perfect and some events are under-or overestimated. However, a majority of the events, including events from the minority classes (impact and overwash) are well classified. The model is then used to predict storm impact hourly on the period without observations . For the period 1993-2016, this model has predicted 232547 Swash events, 741 Impact events and 260 Overwash events.

Bayesian network 4.2.1. Structure learning

Figure 3 shows the BN structures for the two datasets found by the Bayesian search implemented by GeNIe Modeler. This score-based method aim to find the structure that maximizes the out-of-sample classification accuracy while taking into account a list of forbidden arcs. For this study, we forbid all the arcs that are not logical in a physical point of view: storm impact regime can not influence any of the variables, atmospheric surge can not influence any of the variables expect storm impact regime, the waves parameters can not influence the weather conditions and finally the astronomical tide can not be influenced by any of the variables.

Both BNs are similar in terms of structure. Indeed, the storm impact variable is influenced by the tide and the significant wave height in both networks. This is coherent with the literature as the tide and waves are known to play a significant role in coastal flooding (Yang and Liu, 2020). In the network (b), two other variables seem to influence the storm impact regime: wind and wave direction. About the contribution of these two variables, no references have been found, the dependencies between these variable and the storm impact must have been kept for a question of prediction accuracy. Concerning atmospheric surge, both BN found atmospheric pressure and significant wave height as influencing variables. This is also coherent with the literature as the storm surge is the results of the interaction between atmospheric pressure (Harris, 1963), wave characteristics (Bertin et al., 2015) and wind conditions (Arnaud and Bertin, 2014).

Evaluation of the performances

Once the structure and the parameters of the BNs are learned, we can perform the predictions. To that end, the networks are conditioned with the observations of explanatory variables and they output the probabilities of belonging to the different storm impact regimes. For a set of observations, the predicted storm impact regime is the one with the highest probability of belonging. Table 5 presents the descriptive and predictive performances of the BNs. In both cases, the second BN, built on observational and predicted data, shows better performances than the first BN built only on observational data. In the descriptive case, the second BN has better precision and recall resulting in a F 1 -score 21.5% higher than the first BN. The confusion matrices of both networks are presented in appendix ( Table A6). Concerning the predictive performance, the precision of both BN are nearly similar, however the recall is higher for the second BN which results in a F 1score 10% higher for the second BN. Confusion matrices computed on the test set are presented in Table 6. In general, the minority classes (Impact and Overwash) tend to have higher error rate for both BNs. This is expected as minority classes are more difficult to represent as they contain fewer examples than majority classes. It is worth noting that even if the second BN have better predictive performance it still has false positives and false negatives (underestimation or overestimation) for the Impact and Overwash class. This can be problematic for a potential use in an operational EWS.

Discussion

In this section we reflect on a number of aspects related to the advantages and limitations of (i) the prediction of coastal flooding risk with BN based on observational data and (ii) the extension of the database with SLMs.

The main advantage of training BNs exclusively on observational data is that we avoid potential biases related to data generated by process-based models. Indeed, process-based models are not perfect and bias might come from different sources: simplifying assumptions, errors in input data, discretization of the domain (Babovic et al., 2001(Babovic et al., , 2005)). This is especially true for spectral wave models which are known to underestimate wave parameters during storm events (Rakha et al., 2007;Moeini et al., 2012;Arnoux et al., 2018). However by using exclusively observational data, we depend on the precision and data availability of the diverse monitoring systems.

The BN based only on observational data is able to predict coastal flooding risk in a qualitative manner. However, the predictive accuracy of this BN is questionable especially for the minority regimes (impact and overwash) which are the most interesting regimes in terms of coastal flooding prediction. This result was expected as we try to predict extreme events with a database covering only 4 years. For an implementation in a EWS, the predictive accuracy of the BN should be further investigated. More attention could be given to probabilities of belonging for each regime instead of studying only the most likely regime. In addition, sensitivity analysis on decision thresholds could be performed to minimize the number of false positives and negatives for impact and overwash regimes. This step is essential for operational implementation as false positives could lead to the unnecessary implementation of costly mitigation measures such as the installation of sandbag protection while false negatives could lead to significant damages with the absence of mitigation measures.

Concerning the extension of database with SLMs, it led to better results in terms of metrics, nevertheless the BN trained with the extended database still presents false positives and negatives for overwash regime which can be problematic for the use of this BN in an operational EWS. Here again, sensitivity analysis on decision thresholds must be done to minimize the number false positives and negatives. The only pitfall of the extension of the database by SLMs concerns the storm impact model. Even though this model has been chosen due to its generalization performance, it has not been validated yet on historical events. Therefore, we do not know the performance of this model in the historical reconstitution of the storm impact variable. The validation of this model could be performed with databases made from press or insurance archives such as the database proposed by Abadie et al. (2018).

The main limitation of the proposed methodology is related to the storm impact data which are too limited. The proposed extension of the database by SLMs only allows a better representation of the historical events. However, a good predictive model should be able to predict unseen storm events. The observations should be completed by synthetic events and their modeled impact by processbased in the same spirit as Poelhekke et al. (2016). Databases from archives and insurance could be also integrated in the bayesian networks (Abadie et al., 2018) to complete the storm impact data. Another limitation is the inherent assumption of time invariance of the coastal flooding process. Indeed, by not taking into account the change in the site morphology, we consider that the conditions leading to coastal flooding do not differ over time. A last limitation lies in the fact that this model do not predict the coastal flooding risk on the entirety of the Grande Plage of Biarritz but only on the transect used to create the timestacks from which the qualitative storm impact data is extracted.

Conclusion

The aim of this study was to assess the performance in coastal flooding prediction of a BN exclusively based on data measured by diverse monitoring net-works. Because observational data about storm impact regime and atmospheric surge were limited, we proposed a methodology based on statistical learning methods to extend the data about these variables. This methodology was based on cross validation and aimed to select the statistical learning method and associated hyperparameters with the best generalizing ability. Two BNs were trained, one exclusively on the observational data and one with both observational and predicted data. To compare the two networks, their performances were evaluated on the same events.

We demonstrated that it is possible to predict coastal flooding risk in a qualitative manner with a BN based only on observational data with a F 1 -score of 0.628. However, the predictive skill of this network is questionable for the most intense storm impact regimes which are impact and overwash regime. Storm impact and atmospheric surge data were both extended by random forest method which is the method that showed the best generalizing ability in the two cross validation. The extension of the database led to a better BN in term of predictive skill, with precision, recall and F 1 -score in average 20% higher than the BN trained only on observational data.

Even though the predictive skill of the two BNs on the most likely class is not perfect, they can be without a doubt great tools in the prediction of coastal risk. For an operational implementation, it would be more interesting to look at the probabilities of belonging instead of looking at the most likely regime. In addition, a sensitivity analysis could be performed on these probabilities to select tailored decision thresholds for the study site. In any case, the observational data should be completed with synthetic events and their impact simulated with process-based model in order to predict potential unseen events.

Conclusion

In this chapter, we showed that BN can predict coastal flooding in a qualitative manner using exclusively observational data collected by diverse monitoring networks. We also demonstrated that the extension of the database by SLM leads to better predictive performance at the expense of more false positives for "overwash" regime.

Before implementing the presented BN in an operational EWS, further investigations must be made. A sensitivity analysis on the probabilities of belonging must be performed to select tailored decision thresholds for the study site. Concerning the extension of the database with SLM, the performance of the storm impact model must be investigated. The presented storm impact models could be extended to incorporate the temporal dynamic of coastal flooding by using Dynamic bayesian networks (DBN) which are BN including the concept of time. Finally, an interesting work could be to compare the performances of a storm impact model based on BN or DBN with the performances of a storm impact model based on process-based models.

Conclusion

This thesis proposed innovative methodologies based on SLM which contribute to the improvement of coastal risk assessment tools and to the development of an early warning system which aims to reduce coastal flooding risk.

Firstly, we employed SLM to model two coastal processes involved in coastal flooding namely: atmospheric surge and wave runup. We showed that supervised models can be used to predict accurately the coastal processes (in the training data variability) and that knowledge can be acquired by performing variable importance analysis on these models. We also highlighted the ability of unsupervised learning methods to detect groups and patterns from big data, which helped in the characterization of the local sea state and its seasonality. Predictive models and knowledge about coastal processes at a local scale are essential in the development of EWS and consequently in the DRR process.

Secondly, we proposed a methodology based on SLM to improve the forecast of process-based models that are commonly used in operational settings.

We focused more particularly on the improvement of the wave forecast made by spectral wave model which are known to underestimate wave parameters during storm events. With local data and a data assimilation method based on SLM, we were able to improve significantly the forecast of wave parameters during stormy conditions. By improving the wave forecast, this methodology contributes to the improvement of the EWS predictions. This improvement is all the more important because the risk of coastal flooding is at its highest during storm events.

Thirdly, we proposed a new methodology based on deep learning methods to collect storm impact data routinely. We employed Convolutional neural networks (CNN) to classify images from the video monitoring station into storm impact regimes which are categories of coastal flooding risk. Once the CNN are trained, they can be employed to classify the newly created images and therefore generate incrementally a storm impact database. Storm impact databases are crucial in the DRR process as they are used to develop and validate storm impact model or EWS.

Finally, we developed a predictive model translating offshore hydraulic boundary conditions into onshore hazards with SLM. This storm impact model, based on BN, was trained using only observational data collected by different monitoring networks and was able to predict qualitatively the storm impact. In addition, we showed that SLM can be employed to extend the database which can results in a BN with better predictive skill. This storm impact model can be included in an EWS to complement current storm impact model. It has the advantages to avoid the computation time and the potential bias related to the process-based models commonly used.

Scientific contributions

The objective of this research work was to propose innovative solutions to the different scientific issues discussed in each chapter.

The first chapter acted as an "introduction" about the use of SLM in the study of coastal processes. The aim of this chapter was to highlight the ability of SLM to provide better knowledge and understanding about coastal processes at a local scale. Our contribution is related to the fact that we employed methodologies that have never been applied to our study site. In the previous applications of SLM for the modeling of coastal processes, the prediction accuracy was the main focus and little attention was given on the explicability of the SLM. This is why, we also decided to deal with the explicability of the models by performing variable importance analysis.

The second chapter focused on the underestimation of wave parameters by spectral wave model during storm events. This problem is well known in literature and many solutions have been already proposed, including the data assimilation method we employed. However, shallow neural networks were the off-the-shelves SLM in all the previous works using the error prediction method. In addition, little attention was given to hyperparameter tuning in these works. Consequently, we proposed to use ensemble methods (random forest and gradient boosting trees) as alternatives to the neural networks and we performed hyperparameter tuning with Bayesian optimization. With this work, we proved the importance of comparing different statistical learning methods and choosing the optimal hyperparameters in order to obtain the best results.

The third chapter concerned the generation of storm impact databases.

These databases are essential in the DRR process because they are used to train storm impact models or to validate EWS. However, these databases are rare, sparse and mostly come from archives or insurances data. Therefore, we proposed a methodology to collect routinely data about storm impact. This innovative methodology is based on convolutional neural networks (CNN) and images created by video monitoring stations. It aims to classify the images into different storm impact regimes which are categories of coastal flooding risk. To find the best practices for this classification task, we compared several CNN architectures and methods to cope with class imbalance which is a problem related to the extreme nature of storms. We also investigated on the transferability of this method by looking at "pretraining". We showed that a CNN pretrained on a study site can lead to better classification results in fewer epochs for a new study site. Therefore, we shared our pretrained CNN and code with the community to facilitate the application of this methodology to new sites.

The last chapter was about the development of a storm impact model with bayesian network based exclusively on data acquired with diverse monitoring networks. Even though BNs have been employed extensively in the development of storm impact models, they all rely on process-based models or empirical formula to generate the training data. This can be problematic as process-based models and empirical formula are not perfect and can be biased. This chapter aimed to address this problem by building a storm impact model for coastal flooding prediction with a BN exclusively on observational data.

In addition, we proposed and tested a methodology based on SLM to extend the storm impact data as they were limited for our study site.

Finally, a great emphasis was placed on the dissemination and reproducibility of our scientific work. We developed functions in the R package rlm-DataDriven related to an article published at the beginning of this thesis (Appendix B). In this article, we present a robust estimation procedure for auto-regressive models with heterogeneity (Callens et al., 2020). Concerning the improvement of the wave forecast, the code and a short tutorial were made available on a GitHub repository (Link to the tutorial). For the generation of the storm impact database, the code, the weights of the CNN and even the code of the shiny application developed to annotate the timestacks were made available on Github repositories (Link for code and weights of CNN, Link for the application).

Perspectives

This research work opens up a number of perspectives that we have mentioned on different occasions in this manuscript. We extend on some of these perspectives below.

Extreme nature of storm

The main limitation encountered in this thesis is related to the extreme nature of storms. Extreme events are rare by definition and even if the monitoring networks have been recording for years or decades, we always have a limited number of observations of these events. This is a problem in the training of SLM as the storm events, the events of interest, are under-represented in the training dataset. Methods to deal with class imbalance were employed

Operational implementation of the methods

So far, none of the proposed methods in this thesis has yet been implemented in an operational context. Therefore, the operational implementation of these methods is the logical continuation of this thesis. During this step, numerous interrogations will emerge. The first interrogation is about the frequency of "retraining" of the SLM. Because the monitoring networks keep collecting new data, it is legitimate to know the frequency at which the SLM need to be "retrained" to include the information of the newly collected data. This frequency will surely be different for each method proposed in this thesis. Concerning the generation of the storm impact database with CNN, the new timestacks will need to be labeled in order to "retrain" the CNN. To facilitate and even avoid the tedious labeling step of the new timestacks, a special attention could be given to semi-supervised learning which aim to jointly learn from labeled and unlabeled samples [START_REF] Baur | Semi-supervised deep learning for fully convolutional networks[END_REF]. Concerning the storm impact model, a sensitivity analysis on the probabilities of belonging must be performed to select optimal decision thresholds. This sensitivity analysis is crucial for operational implementation as it aims to minimize the false positives for the storm impact regimes representing the highest risk.

A. Extracting storm surge data with harmonic analysis

The observed storm surge is obtained by subtracting the astronomical tide from the measured water level of a tide gauge. The astronomical tide has been studied for a long time and many methods allow for the prediction of the tide level. Harmonical analysis is the most commonly used method and the one presented below. This method aims to represent the tidal water level as a sum of basic harmonic constituents.

Data

The water level data from the tide gauge located in Socoa (city in the south west of France) were provided by the french Naval Hydrographic and Oceanographic Service (SHOM). The data range from 2011 to nowadays with a hourly time step. In addition to provide the data, the SHOM also documented each dysfunction of the tide gauge. This allowed for the proper removal of aberrant data. The final data are presented in the

Method

Harmonic analysis is the most common method to model astronomical tide with water level from tide gauge. In this method, the tidal level is modeled as a sum of harmonical constituents, the majority of which being related to the gravitational forces of the Moon and Sun [START_REF] Nicolle | Modélisation Des Marées et Des Surcotes Dans Les Pertuis Charentais[END_REF]. At time t, the tide level can be written as :

h(t) = Z 0 + n i=1 [A i cos(v i t + φ i )], (A.1)
where Z 0 corresponds to the mean water level, A i , v i , φ i to the amplitude, the frequency and the phase of the wave respectively. Finally, n corresponds to the number of constituents chosen.

The unknown parameters A i and φ i can be estimated with the least square method for each constituent. An implementation of this method in R code is available in the package TideHarmonics. For the harmonic analysis method, the number of harmonic constituents must be chosen: a higher number means a better accuracy. Because we want the most accurate model, we chose to estimate the coefficient of the first 114 harmonic constituents which is the For the modelling of the tidal water level, there will always be some discrepancies between the predicted and the measured values at the tide gauge.

Indeed, the tide gauge measures the total water level which is composed by the addition of the astronomical tidal level and the water level which is induced by the meteorological conditions. The low pressure and high winds during a storm lead to an increase in water level which is called "storm surge".

Observed storm surge can be extracted by subtracting the modelled astro-

Conclusion

Harmonical analysis was used to model the tide with the water level data from the tide gauge data of Socoa. It was found that the type of the tide in this area was semi-diurnal. In addition, this method allowed us to compute the observed storm surge for the site of Socoa by subtracting the modelled astronomical tide from the water level observed at the tide gauge. In environmental modelling, pure homoscedasticity is uncommon. For example, residuals of hydrological [8], or air pollution models [26] are usually heteroscedastic. Ignoring this problem and performing an ordinary least squares regression would result in regression parameters with biased covariance matrix and hence would lead to erroneous inference. One method to deal with this biased covariance matrix is to use the White's estimator [25] which provides a heteroscedasticity consistent covariance matrix. In case of a time series, the Newey-West estimator [17] can provide heteroscedasticity and an autocorrelation-consistent covariance matrix. Both of these methods account for heteroscedasticity but do not give information on the variability of the data generation process.

Another approach to cope with heteroscedasticity is to perform a weighted analysis, where the assumed underlying model is:

y i = x T i β + σ i i ,
with y i the observed response, x i the associated covariates, β the vector containing the parameters to be estimated, i the independently and identically distributed error terms with mean 0 and unknown symmetric distribution function and σ i the term accounting for heteroscedasticity. Usually, a parametric function is assumed for this term [1,4,7,9]. It can be a power function of the mean as proposed by Box et al. [3]: σ i = φ|x T i β| γ or some functions of the covariates. Parameters in the variance function are not known and have to be estimated by maximum likelihood method.

Estimation of parameters in models exhibiting heterogeneous variance is performed by an iterative procedure. A preliminary estimate of mean parameters is obtained by the least squares method. Residuals of this model are then used to estimate variance parameters. Finally, a weighted least squares method is performed with the estimated variance as weight. Unlike White and Newey-West methods, modeling heterogeneous variance enables one to get better estimates for the mean parameters and also to gain information on the variability of the data generation process [7,9]. Such as least squares method, the estimation method for heteroscedastic models has a low breakdown point of 1/n, meaning that only one outlier in the observed response can have a large effect on the estimation of the mean parameters [20,23,27]. The estimation of the variance parameter is also affected as maximum likelihood methods are very sensitive to outliers [9,21].

In practice, outliers are prevalent in environmental dataset. They can be found in both response variable and covariates but, in this work, our interest lies only in outliers that are present in the response variable. In the presence of outliers, robust methods must be applied. They aim to produce reliable estimates that are not seriously affected by outliers, extreme values or small deviations from model assumptions [13].

When both heteroscedasticity and outliers are present in regression analysis, one can use the method described by Carroll et al. [4]. This iterative method allows one to robustly estimate mean and variance parameters. For the mean parameters, they perform a weighted M-estimation with variance as weight. M-estimation is a robust method which consists in minimizing a loss function that is slowly varying for abnormal residuals instead of squared residuals [24]. This loss function is controlled by a tuning parameter c which "regulates the amount of robustness" [11]. For the variance, they assume a parametric function and robustly estimate its parameters with high-breakdown point estimators. A limitation of this method is that errors must be independently and identically distributed, which is not the case when regression is performed on temporal data.

In this article, we propose an estimation procedure which adapts the weighted M-estimation method of Carroll et al. [4] to time series by taking into account temporal correlations, and, chooses the value of the tuning constant by minimizing the variance of the estimators such as the work of Wang et al. [23,24]. We perform numerical studies to compare the proposed method with other usual regression methods. Finally, we illustrate our methodology with a dataset on fine particle matter pollution (PM 2.5 ) in the city of Beijing. This dataset is particularly interesting for two reasons. First, variability of PM 2.5 concentration is not likely to be constant. Second, largest concentrations, which are related to increased coal consumption during cold days, may have a significant impact on the estimation of regression parameters.

Method

The regression model

Let y = (y 1 , y 2 , ... , y n ) T be the observed response measured over n equivalent time periods and x i = (x i1 , ... , x ik ) T be the set of k associated predictors (x i1 is equal to 1 if intercept is considered as a predictor). We assume data are generated from the following heteroscedastic linear model :

y i = x T i β + σ i i , (1)
in which β is the vector collecting the parameters to be estimated, i are the error terms following an autoregressive process of order p (AR(p)) which takes into account temporal correlations, and

σ i = φg(x T i β, γ),
where g(•) is a known function of the mean (x T β) with unknown parameter vector γ, and unknown dispersion parameter φ. Several choices for σ i are possible, few examples are : * σ i = φ(1 + |x T i β|) γ or σ i = φ|x T i β| γ proposed by Box et al. [3], * σ i = φe γx T i β proposed by Bickel et al. [2].

Estimation of the parameters

For given or estimated value σi , the robust M-estimation for β minimizes :

n i=1 ρ y i -x T i β σi ,
where ρ is a loss function that is slowly varying for abnormal residuals (outliers). The most known is Huber's loss function :

ρ(u) = 1 2 u 2 if |u| ≤ c c|u| -1 2 c 2 if |u| > c . (2)
Here, c is a tuning parameter chosen between 0 and 3, which controls the degree of robustness. Default value for Huber's function is 1.345 to ensure 95% asymptotic relative efficiency when data are normally distributed. More examples of loss functions can be found in Wang et al. [23].

Taking derivatives of (2) leads to the following estimating equation of β:

U (β) = n i=1 x i σi ψ y i -x T i β σi = 0,
where ψ(x) = min{c, max{x, -c}} is the derivative of Huber's loss function.

To solve this estimating function, one can rewrite U (β) as :

U (β) = n i=1 x i W i r i = 0,
where W i = ψ(r i )/r i are weighting terms, and r i = (y ix T i β)/σ i are the Pearson residuals. Now, for a given weight W i , the robust estimator of β can be obtained by the following formula :

β = { n i=1 x i W i x T i } -1 { n i=1 x i W i y i }. (3)
An iterative approach is needed as W i is a function of β and σ. This approach is derived from the pseudolikelihood approach and consists in fixing alternatively parameters of the variance (γ and φ) and the regression parameters (β).

The variance parameters are also robustly estimated and are given by :

-A high breakdown estimator for γ :

n i=1 χ y i -x T i β φg(x T i β, γ) g (x T i β, γ) g(x T i β, γ) = 0, (4)
where χ(•) is a bounded function. Croux [6] and Bianco et al. [1] used χ(x) = min(x 2 /c 2 1 , 1) -0.5 with c 1 = 1.041 to obtain a 50% breakdown estimator of γ under the normality assumption.

-The MAD estimator for the dispersion parameter : φ = Median |y ix T i β| g(x T i β, γ) /0.6745.

Wang et al. [24] showed that under some regularity conditions, the robust estimator β obtained by the iterative procedure is Fisher consistent. Moreover, when n → ∞ the covariance matrix is given by [4,10,11] :

var( β) = K 2 [1/(n -k)] n i=1 ψ(r i ) 2 [(1/n) n i=1 ψ (r i )] 2 × n i=1 {σ 2 i (x T i x i ) -1 },
where r i = (y ix T i β)/σ i , and

K = 1 + p n var(ψ (r i )) (Eψ (r i )) 2 .
with p the number of unknown parameters.

A data-dependent tuning constant

The tuning parameter c associated with the loss function regulates the amount of robustness in the estimation. When observations are normally distributed and without outliers, optimal value of this parameter should be c = +∞. On the other hand, for heavy tailed distributions, optimal value should be around or smaller than 1. The value of the tuning parameter should be chosen carefully since robustness comes at the price of efficiency. Indeed, a smaller value than needed, would result in considering more usual observations as outliers and would lead to a loss of efficiency in the estimation of the regression parameters.

As in Wang et al. [23], we define the best tuning constant as the one which minimizes the variance of the regression parameters. Therefore we propose to repeat the estimating procedure with different values of c between 0 and 3 for the Huber's function [23] and choose the one which minimizes the sum of the estimated variances of the regression parameters.

Accounting for temporal correlations

So far we have only considered the independent model, we now need to incorporate the autoregressive process of order p present in the error terms. We write i as p j=1 (α j i-j ) + ξ i where ξ i are independent errors and we rewrite the model (1) as :

y i = x T i β + p j=1
α j σ i i-j + σ i ξ i .

Because the i are unobserved, we propose to use the Pearson residuals from the initial model ( 1), say, εi , and we now have the following linear model with roughly independent errors,

y i = x T i β + p j=1
α j σi εi-j + η i , where ( σi εi-1 , σi εi-2 , ..., σi εi-j ,) are the augmented additional covariates including p lagged residuals, and (β, α 1 , ..., α p ) are the new parameters to be estimated including p lag parameters, σi is an estimate of σ i . In the iterative procedure to be described below, this σi will be estimated from the variance function using the previous parameter estimates for (φ, γ, β). Here, η i represents the resulting error which should be close to σ i ξ i . We fit this augmented model with the optimal value of the tuning parameter to obtain the final estimate of β.

In the application section, we will demonstrate how we choose the order p of the autoregressive process.

The estimation procedure

The complete estimation procedure is summarized in the following algorithm :

1. Obtain an initial robust estimate β0 assuming a constant variance g(x T β, γ) = 1 and using M-estimation with the default value of c (rlm function).

2. By fixing β = β0 , the robust variance parameters ( φ, γ) are estimated with ( 4) and ( 5) respectively.

3. By fixing the variance parameters equal to their robust estimates, we update β with (3).

4. To find the best tuning parameter, the steps 2-3 are repeated for a range of c values between 0 and 3. The best tuning constant is the one which minimizes the sum of the estimated variance of the regression parameters.

5. The model is fitted using the best value of the tuning constant (ĉ) from the steps 2-3. Then, temporal correlations are added by following the procedure described previously.

Numerical studies

In this section, we investigate the performance of our procedure. We compare mean bias and mean square errors of the estimates obtained by different methods such as least squares (lm function in R), generalized least 4 Application to Air Quality Data

Context of the application

We apply the proposed robust procedure to analyse fine particle matter (PM 2.5 ) concentrations from the US embassy situated in Beijing, China. Such as other China's mega cities, Beijing has been suffering from chronic air pollution [5]. The main constituents of this pollution are suspended particle matters under 2.5 µm of diameter, widely known as PM 2.5 [15]. Concentration of PM 2.5 is highly variable and depends on sources of emission, secondary chemical generation processes and meteorological conditions. According to numerous studies [15,18,28], these suspended particles affect climate, visibility and human health in many ways.

Given the severity of the pollution and the potential hazardous effects, China's State Council aimed to reduce the PM 2.5 pollution by at least 25% for the period 2012-2017. Our objectives are (i) to show the efficiency of our procedure on a real dataset compared to a conventional estimation method (least squares method), (ii) to analyze the potential effect of the decision of China's State Council on the concentration of PM 2.5 . Our method is highly desirable in this case study since the variance of PM 2.5 concentrations is not likely to be homogeneous [26] and outliers may be present.

Regression analysis

The hourly data come from a previous study lead by Liang et al . [15] (Figure 2). PM 2.5 concentrations were taken at the US Embassy of Beijing and meteorological measurements at the Beijing Capital Airport. Both time series covers the period from 1 January 2010 to 31 December 2014.

To evaluate changes in PM 2.5 concentration after the decision of China's State Council, we created two new variables : Policy and Time Policy. Both variables take the value 0 before 2012, however Policy take the value 1 after 2012 to detect any shifts in the intercept and Policy Time take the value of the time lapsed in days after 2012 to test and quantify the trend after that decision.

We modelled PM 2.5 concentration as a linear combination of the following covariates : dew point ( • C), temperature ( • C), atmospheric pression (hPa), combined wind direction (3 factors), cumulated wind speed (m/s), cumulated hour of snow (mm), cumulated hour of rain (mm), Time policy, Policy and seasonal patterns with sin and cos functions for the 3,2 years and 6,4,3 months cycles. To perform the regression analysis, the following numerical covariates have been standardized: dew point, temperature, atmospheric pression, cumulated wind speed, cumulated hour of snow, cumulated hour of rain.

The response variable PM 2.5 is a concentration; therefore a transformation was needed in order to avoid predictions lower than zero with the regression analysis. Two transformations were tested on the dataset: square-root and logarithm transformation (Figure 2). For both transformations, mean squared and mean absolute errors were computed after re-transformating fitted values in original scale. Althougth logarithm transformation is the classical transformation for air quality data, we chose the square root transformation as it yielded slightly lower errors (Supporting information, Table S-3).

The model is fitted with least squares method by the lm function from R statistical software [19]. The residuals vs. fitted value plot (Figure 3) indicates the presence of heteroscedasticity with larger residuals as fitted values increase. This heteroscedasticity does not lead to biased estimators but to estimators with biased covariance matrix. This could result in underestimation of standard errors, erroneous Z-values and therefore erroneous hypothesis tests. The normal probability plot in Figure 3 reveals that distribution of residuals is skewed, indicating the presence of outliers which may have influenced the estimation of regression parameters. Moreover, temporal correlations in the residuals have been found using the ACF and pACF plots.

In this present case, our method is highly desirable as we have the presence of heteroscedasticity, temporal correlations and outliers.

Robust regression analysis with the proposed method

We fitted the same regression model to the data with the proposed method. This method has been implemented in a R function: rlmDD het, available in the package RlmDataDriven. The R code for the presented analysis can be found in the supporting information file.

In the literature, we did not find previous papers or indications on how to model variance of PM 2.5 concentration, consequently, we used common variance functions for the analysis such as power or exponential functions. Hereafter, we present the result obtained by considering σ i = φe γ|x T β| .

First, our method chooses best tuning constant by testing a range of values for the tuning parameter between 0 and 3. As stated earlier, the best value of c will be the value that minimizes the sum of the variance of the regression parameters. For the PM 2.5 data, optimal value of the tuning parameter is found around 1.5.

Then, we use ACF and pACF plots of the robustified residuals to determine the order of the autoregressive process. The robustified residuals are defined as ψ((y ix T i β)/σ i ). Two significant lags are found in the pACF plot, therefore, we consider that p = 2 and we add two lagged terms to the regression model. The lagged terms are: α 1 σi εi-1 + α 2 σi εi-2 + η, where σi εi-1 , σi εi-2 are lagged terms built from the initial model estimated with the best tuning constant (ĉ = 1.5) and η is the assumed independent error term. The term ε corresponds to Pearson residuals of the initial model and σ is the estimated variance function.

The normal probability plot (Figure 3) of the residuals clearly illustrates that our robust procedure has taken care of outliers successfully. The residuals vs. fitted values plot (Figure 3) seems to indicate that heteroscedasticity has been lowered.

The results of both methods are listed in Table 2. In this table, the covariance matrix of the regression parameters obtained by least squares method has been estimated with the NeweyWest function (sandwich package) which gives a heteroscedasticity and autocorrelation consistent estimation of the covariance. This estimation was necessary to obtain corrected standard errors as the residuals exhibited heteroscedasticity and temporal correlations.

Our estimation method drastically reduced the variance of the parameters. It is worth noting that the covariates Policy and Time policy are not significantly different from 0 in the least squares method contrary to our robust method. Years after 2011 are characterized by an positive shift in the intercept. However, coefficient of the variable Time policy is slightly negative. After one year, this coefficient outweights the positive shift in the intercept, meaning that the PM 2.5 concentration decreases slightly compared to years before 2012.

Finally, we can see in the Table 2 that lagged terms are significant. This indicates that the two previous terms contribute significantly to the output and were, therefore, necessary to consider.

Discussion

This method is data-dependent by the optimal choice of the tuning constant and it incorporates temporal correlations by adding lagged terms in the covariates. The numerical study showed that this procedure outperforms the other usual regression methods when data are contaminated by providing more precise estimates for the mean parameters. In the application with the PM 2.5 concentration dataset, we proved that our method results in estimates with significantly lower variance compared to the ones obtained by least squares estimation, leading to better hypothesis testing. This method is well suited for environmental dataset due to the frequent presence of heterogeneity and outliers. Hyslop et al. [12] utilized Thiel-Sen robust regression to evaluate long-term trends in aerosol concentrations via the historical PM2.5 element measurements. Van Donkelaar et al. [22] used a geographically weighted regression (GWR) statistical model to represent bias of fine PM2.5 concentrations over North America . Knibbs et al. [14] utilized the land-use regression (LUR) to estimate PM2.5 at continental scale and explained the most spatial variability in PM2.5 in Australia . These three methods did not consider the heterogeneity and autoregressive errors. In the future, the proposed method could be generalized to these three models. Furthermore, the proposed method could be extended to model time-series of counts by using a link function such as the generalized linear model [16]. Concerning the variance, it would be helpful to consider more flexible approaches for the estimation of the parameters and to provide guidance for choosing the most appropriate variance function for a dataset.

Supplementary Material

Online supplementary material includes :

-Table S-1 : Mean bias from true parameters and standard errors of the estimators obtained by several regression methods for the case of the power variance function.

-Table S-2 : Number of non convergence removed for the calculation of the mean bias and standard deviation for the simulation study.

-Table S-3 : Mean Absolute Error and Root mean squared error for different transformations of PM 2.5 and different models.

-Figure S-1 : Mean square errors of the estimated parameters obtained by several regression methods for the power variance function.

-Page 5 : Information on the dataset with downloading link and indications on the different steps to tidy the dataset before analyzing it.

-The R script used to obtain the results in the application section.

2. 3 . 2

 32 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Wave Climate Characterization with statistical downscaling . 2.4.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Statistical downscaling . . . . . . . . . . . . . . . . . . 2.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Improve wave forecast at a specific location with ensemble methods and local observations 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Article: Using Random forest and Gradient boosting trees to improve wave forecast at a specific location . . . . . . . . . . . 3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Automatic creation of storm impact database based on video monitoring and convolutional neural networks 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Article: Automatic creation of storm impact database based on video monitoring and convolutional neural networks . . . . . . 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Bayesian networks to model storm impact using data from both monitoring networks and statistical learning methods 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Article: Bayesian networks to model storm impact using data from both monitoring networks and statistical learning methods 5.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Map representing the population density for Europe in 2018. This map was reprinted from ec.europa.eu/eurostat.

Figure 1 . 2 :

 12 Figure 1.2: The concept of risk.

Figure 1 . 3 :

 13 Figure 1.3: Coastal flooding on the Grande Plage de Biarritz during the Christine storm (04/03/2014). Image reproduced from www.rtl.fr.

Figure 1 . 4 :

 14 Figure 1.4: Different components of total water level. Reprinted from Vitousek et al. (2017).

Figure 1 . 5 :

 15 Figure 1.5: Disaster reduction cycle. Reprinted fromYang and Liu (2020)

Figure 1 Figure 1 . 6 :

 116 Figure 1.6: Normalized measure of the interest made by Google Trends for different terms related to statistical learning. Measures are represented in grey lines and moving average over 6 months are represented in red lines.

20 2

 20 The objective of this chapter is to demonstrate the interest of SLM in the modeling of coastal processes by showing that they can provide knowledge about these processes or the study site. This chapter presents three examples where SLM were employed to model two different processes which are storm surge, wave runup on the beach and to characterize wave climate.In the first part of the chapter, we show how to model storm surge level with several environmental variables related to the tide, the winds and the waves. Different methods are tested to find the best predictive model and variable importance analysis is performed on the best model to learn which variable is the most predictive for the study site. The insights brought by the variable importance analysis are discussed with the characteristics of the site and what is known in the literature.

Figure 2 . 1 :

 21 Figure 2.1: Map showing the locations of the directional wave buoy (1), the meteorological station of Biarritz (2) and the tide gauge of Socoa (3).

Figure 2 . 2 :

 22 Figure 2.2: Atmospheric surge level extracted from the water level of Socoa tide gauge.

Figure 2 . 3 :

 23 Figure 2.3: Permutation importance of variables used in Random Forest to predict storm surge level

Figure 2 . 4 :

 24 Figure 2.4: Measured runup on the beach versus the runup predicted by Stockdon formula and the local formula computed for the GPB.

Figure 2 . 5 :

 25 Figure 2.5: Transects monitored on the Grande Plage de Biarritz and example of associated timestack. The red line in the timestack represents the position of the seawall.

Figure 2 . 6 :

 26 Figure 2.6: RMSE and associated standard error obtained for the 10 repeated 10-fold cross validation for different statistical learning methods.

Figure 2 . 7 :

 27 Figure 2.7: Permutation importance computed for the random forest model.

Figure 2 .

 2 Figure 2.9 shows the smoothing functions for the GAM model. For the significant wave height H s the smoothing function is not monotonic. It increases as the value of H s increases, however above 3.8 meters, the function starts to

Figure 2 . 8 :

 28 Figure 2.8: Basic model checking plots for GAM fitted on the data. The upper left panel represent a QQ-plot of the residuals while the upper right panel represent the residuals vs the linear predictor. The lower left panel is an histogram of the residuals and the lower right is a plot of the fitted values against the residuals.

Figure 2 . 9 :

 29 Figure 2.9: Spline functions estimated by the GAM model.

Figure 2 . 10 :

 210 Figure 2.10: Permutation importance computed for the GAM model.

  predictions of the models (fitted on the whole dataset of the measurement campaign) plotted against H s are presented in Figure2.11.

Figure 2 . 11 :

 211 Figure 2.11: Plot of the runup predictions made with the buoy data against H s . The red dots are the data collected during the measurement campaign and served as training data.

  rection at peak θ p are recorded by a directional wave rider buoy (DWR MKIII) operated by the Center for Studies and Expertise on Risks, Environment, Mobility, and Urban and Country Planning (CEREMA) and the University of

Figure 2 . 12 : 1 •

 2121 Figure 2.12: Spatial distribution of seal level pressure data with resolution of 1 •

  step. These 26 components explains 99% of the total variance of the original data. The data projected on the first 26 components of the PCA are mapped on a 10 × 10 rectangular grid by the SOM algorithm. A sensitivity analysis based on the total within-cluster sum of square (wss) was performed to select the number of cluster. The number of 100 clusters (10 × 10) was chosen as it showed great performance in term of wss and as this number was large enough to capture a wide range of atmospheric situations.The mean SLP fields of the 100 weather types obtained by the SOM algorithm are displayed in Figure2.13. From this figure, we can see that similar weather types are closer on the bidimensional lattice. On the upper-left corner, the majority of weather types display high-pressure systems in the Atlantic Ocean. On the contrary, on the lower-left corner the weather types mostly display intense low-pressure systems. These weather types can be identified as the two mains circulation pattern in the North Atlantic Ocean: North Atlantic oscillation (NAO) and the east Atlantic oscillation (EA). The former circulation pattern (NAO) is characterized by the presence of both a low-pressure system located over Iceland and a high-pressure pattern located around the Azores area[START_REF] Mellado-Cano | Examining the north atlantic oscillation, east atlantic pattern, and jet variability since 1685[END_REF]. This circulation pattern is represented in the upper part of the lower-left corner in the lattice. The latter circulation pattern (EA) is characterized by the presence of a low-pressure system in the south of Iceland and west of Ireland[START_REF] Mellado-Cano | Examining the north atlantic oscillation, east atlantic pattern, and jet variability since 1685[END_REF]. The weather types presenting this pattern are located in the lower part of the lower-left corner of the lattice. On the right side, the weather types in the middle display high-pressure system while weathers type on the upper and lower corner show low-pressure systems (not as intense as the left side).

Figure 2 . 13 :

 213 Figure 2.13: Weather types obtained by the self organizing map method. Blue color represents zones of low pressure and red color zone of high pressure.

Figure 2 . 14 :

 214 Figure 2.14: Observed frequencies of the different weather types for the period 1969-2019

Figure 2 . 15 :

 215 Figure 2.15: Observed frequencies of the different weather types depending on the seasons.

Figure 2 .

 2 16 shows the bivariate graphs of H s and T p for each WT. The distribution of H s and T p change significantly between WT. The highest values for both parameters are commonly found in WT displaying low-pressure system whereas the lowest values are found in system presenting strong high-pressure system or homogeneous pressure field. Unlike the significant wave height, the wave direction (θ) do not vary between weather types. For all the observations, this variable do not vary a lot: 75% of the observations have a direction contained in the 52 following interval: [296.6, 312.6]. The Figure2.17 shows the histogram of θ and bivariate graphs of H s depending on θ for the most represented weather types (WT28 and 56).

Figure 2 .Figure 2 . 16 :Figure 2 . 18 :Figure 2 .

 22162182 Figure 2.17: (a) Histogram of the wave direction for all the observations and (b) Bivariate graphs of significant wave height (Hs) and wave direction (θ) for the most common weather types.

59 3 .

 3 Improve wave forecast at a specific location with ensemble methods and local observations 3.1 Introduction

Figure 1 :

 1 Figure 1: Map showing the location of the study site. The red dots show of the locations of the directional wave buoy (1), the meteorological station (2) and the beach called "Grande Plage de Biarritz" (3).

Figure 2 :

 2 Figure2: Distribution of the errors computed between values observed at the buoy and values corrected or not with the different machine learning algorithms. "Num" stands for numerical model (no correction), "Ann" for artificial neural networks, "Rf" for random forest and "Gb" for gradient boosting trees. The horizontal lines in the red boxplots represent from top to bottom: the third quartile, the median and the first quartile.

Figure 3 :

 3 Figure 3: Variable importance for the correction of the three wave integrated parameters.

Aurélien

  Callens a, * , Denis Morichon b , Stéphane Abadie b , Matthias Delpey c , Benoit Liquet a,d a Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France b Université de Pau et des Pays de l'Adour, E2S UPPA, SIAME, Anglet, France c Centre Rivages Pro Tech SUEZ EAU FRANCE, Bidart, France d Department of Mathematics and Statistics, Macquarie University, Sydney, Australia Distribution of θ p measured at the buoy.

Figure 1 :

 1 Figure 1: Wave direction vs significant wave height. Both parameters are measured at the buoy of Anglet.

Figure 3 :

 3 Figure 3: H s values observed at the buoy vs H s values corrected with different machine learning algorithms

Figure 5 :

 5 Figure 5: θ p values observed at the buoy vs θ p values corrected with different machine learning algorithms

•

  [START_REF] Callens | Automatic creation of storm impact database based on video monitoring and convolutional neural networks[END_REF].Automatic Creation of Storm Impact Database Based on Video Monitoring and Convolutional Neural Networks. Remote Sensing, 13(10), 1933. Communication: • "Automatic Creation of Storm Impact Database Based on Video Monitoring and Convolutional Neural Networks", Online seminar for the Probability and Statistics research team of the LMAP (July 2021).

4. 2

 2 Article: Automatic creation of storm impact database based on video monitoring and convolutional neural networks

1

 1

Figure 1 .

 1 The three categories of storm impact regimes estimates from timestack images. y-axis represents the time in seconds, x-axis represent the pixel index (cropped images), and red lines represent the sea wall bottom and top positions. (a) Swash regime. (b) Collision regime. (c) Overwash regime.

Figure 2 .

 2 Figure 2. Map showing the locations of the two study sites.

2. 1 .

 1 Grande Plage de Biarritz 2.1.1. Site Characteristics

Figure 3 .

 3 (a) Satellite view of the site of Biarritz from Google Earth with red lines representing the transects on the site. (b) Transect on the Grande Plage de Biarritz.

Figure 4 .

 4 Figure 4. (a) Satellite view of the site of Zarautz from Google Earth, with the red point representing the location of the video monitoring station. (b) Positions of the transects on the site of Zarautz.

Figure 5 .

 5 Figure5. Workflow for this study. Items inside the boxes that are highlighted in red represent the choices tested in this study, whereas the items in black are the methods applied in every cases.

F 1 =

 1 Precision × Recall Precision + Recall where Precision = True positives True positives + False positives and Recall = True positives True positives + False negatives .

Figure 6 .

 6 Figure 6. Prediction errors on the test set of Biarritz by the best CNN. The storm impact regime predicted is written under each timestack, and the ground truth is written in parentheses. The probabilities of belonging predicted by the CNN are represented on the right side of each timestack (red = prediction made by CNN, blue = ground truth). The red line in the timestack represents the position of the seawall.

  113

  In this chapter, we develop a storm impact model for the Grande Plage de Biarritz based on Bayesian networks (BN). In the literature, most of the applications of BN in storm impact modeling are relying on process-based modeling to constitute the training database for the BN. Unlike previous works on this subject, we only use observational data collected by diverse monitoring networks (tide gauge, wave buoy, weather station) as training database. Because observations of storm impact and atmospheric surge are limited, we also present and test a methodology based on SLM to extend the dataset. This methodology is based on cross validation and aims to select the SLM with the best generalizing ability. In total, two BNs are trained in this chapter, one exclusively on the observational data and one with both observational and predicted data. Their performances are compared by predicting on the same events.

Aurélien

  Callens a , Denis Morichon b , Benoit Liquet a a Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France b Université de Pau et des Pays de l'Adour, E2S UPPA, SIAME, Anglet, France

Figure 1 :

 1 Figure 1: Map showing the location of the study site. The red dots show of the locations of the directional wave buoy (1), the weather station (2), the tide gauge (3) and the beach called "Grande Plage de Biarritz" where a video monitoring station is installed (4).

Figure 2 :

 2 Figure 2: Diagrams of the development of the bayesian networks based on observation data (a) and on observation and predicted data (b).

  T rue positives + False negatives .

Figure 3 :

 3 Figure 3: Bayesian network structures found for (a) the dataset based on observational data and (b) the dataset based on both observational and predicted data.

Figure A. 1 .

 1 To assess the accuracy of the harmonic model on unseen data, the dataset has been divided in two part : the training set corresponding to 80% of the data and ranging from 2011-04-26 to 2016-10-03 and the test set corresponding to the remaining 20% and ranging from 2016-10-04 to 2020-06-12. It is worth noting that the period 2016-2020 contains a lot of missing data due to a dysfunction of the tide gauge (Figure A.1).

Figure A. 1 :

 1 Figure A.1: Tide data from the Socoa tide gauge.

  Figure A.2: Water level measured and modelled for the month of May 2020.

Figure A. 3 :

 3 Figure A.3: Atmospheric surge level extracted from the water level of Socoa tide gauge.

Fig. 1 Fig. 2

 12 Fig.1Mean square errors of the estimated parameters obtained by several regression methods for the exponential variance function (log-scale). Results based on 500 replications. Gls corresponds to generalized least squares, Lm to least squares, RlmDDh to the proposed method, Rlm to M-estimation and Whm to weighted M-estimation.

Fig. 3

 3 Fig.3Regression diagnostic plots for the least squares model (upper panel) and the proposed method (lower panel). In the residuals vs fitted values plot for the method, 6 residuals lie outside the plot area.

Table 2 .

 2 1: Examples of atmospheric surge modeling with SNN in the literature.

	2.2.2 Methodology
	Supervised learning methods
	Over the last few years, the use of Shallow neural networks (SNN) has been
	democratized for the modeling of water level related to the atmospheric surge,

table 2

 2

		.2 shows the RMSE and MAE computed on the test data for the
	different models. From this table we can see that the performances of the
	different learning algorithms are quite close. Random forest shows slightly
	lower RMSE and MAE. The results obtained for our site are consistent with
	the literature. For example, Lee (2006), who performed a similar analysis with
	neural networks on several Taiwanese stations, found RMSE values between
	5.39 and 9 cm.
	Table 2.2: MAE and RMSE (in cm) computed on the test data for the different
	models.
		Linear model Random Forest Gradient boosting trees Neural Networks
	RMSE	6.61	5.73	6.54	6.56
	MAE	4.86	4.09	4.78	4.81

  wave parameters such as significant wave height H s , peak period T p , wave

	directional spreading (disp) and direction at peak θ p are recorded by a direc-
	tional wave rider buoy (DWR MKIII) operated by the Center for Studies and
	Expertise on Risks, Environment, Mobility, and Urban and Country Planning
	(CEREMA) and the University of Pau and Pays de l'Adour (UPPA). The buoy
	is located a few miles off GBP at 50 meters water depth. Since its deployment
	in 2009, this buoy has been recording wave parameters every 30 minutes. By

Table 2 .3: Characteristics of the different statistical learning models

 2

	Method used	Response variable Covariates	Hyperparameters
	Linear models Runup	H s , T p , disp, θ and β	None
	Linear models Runup	Log transformed H s , T p , disp, θ and β	None
				Default smooth functions (thin plates)
	GAM	Runup	Splines for H s , T p , disp, θ but not for β	Additional penalty term in the smoothness selection procedure with select = TRUE in
				mgcv package.
	Random forest Runup	H s , T p , disp, θ and β	Default values (number of trees = 500)

•

  [START_REF] Callens | Using random forest and gradient boosting trees to improve wave forecast at a specific location[END_REF]. Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France b Université de Pau et des Pays de l'Adour, E2S UPPA, SIAME, Anglet, France c Centre Rivages Pro Tech SUEZ EAU FRANCE, Bidart, France d Department of Mathematics and Statistics, Macquarie University, Sydney, Australia

	Using Random forest and Gradient boosting trees to
	improve wave forecast at a specific location
	Aurélien Callens a, * , Denis Morichon b , Stéphane Abadie b , Matthias Delpey c , Benoit Liquet a,d
	Using Random forest and Gradient boosting trees to improve wave fore-
	cast at a specific location. Applied Ocean Research, 104, 102339.
	3.2 Article: Using Random forest and Gradi-
	ent boosting trees to improve wave fore-
	cast at a specific location

Communications:

• "Using Random forest and Gradient boosting trees to improve wave forecast at a specific location", JDS 2021: 52èmes Journées de Statistique de la Société Française de Statistique (SFdS), Online conference (June 2021).

• "Improving numerical wave models with machine learning algorithms", Statistical Modelling in Ecology and Environmental Data Workshop, Anglet (September 2019). a

Table 1 :

 1 Statistical metrics for the three variables of interest before the hyperparameter tuning. "Ann" stands for artificial neural networks, "Rf" for random forest and "Gb" for gradient boosting tree.

	Hs					Tp				θ p
	Numerical	Ann	Rf	Gb	Numerical	Ann	Rf	Gb	Numerical	Ann	Rf	Gb
	model	Corr.	Corr.	Corr.	model	Corr.	Corr.	Corr.	model	Corr.	Corr.	Corr.
						Computed with all data
	Biais	-0.201 0.005 -0.002 -0.004	0.712 0.002 -0.026 0.005	-0.249 0.698 0.976 0.765
	RMSE	0.399 0.306 0.248 0.267	1.839 1.603 1.282 1.388	13.803 12.391 9.749 10.645
	SI	0.166 0.148 0.120 0.129	0.153 0.145 0.116 0.125	0.045 0.041 0.032	0.035
	Cor	0.954 0.962 0.975 0.972	0.78 0.804 0.880 0.857	0.330 0.366 0.455	0.386
	Biais	-0.536 -0.156 -0.124 -0.126	Computed with data where H s > 3m 0.683 -0.026 -0.090 -0.041	1.925 -0.561 0.052	0.046
	RMSE	0.766 0.515 0.420 0.433	1.348 1.083 0.956 1.022	7.351 5.769 4.449 4.931
	SI	0.133 0.120 0.098 0.101	0.089 0.083 0.073 0.078	0.024 0.019 0.015 0.016
	Cor	0.818 0.857 0.902 0.898	0.832 0.850 0.886 0.869	0.589 0.604 0.790 0.726

Table 2 :

 2 Default values, ranges and selected value of hyperparameters for the machine learning algorithms

	Machine learning algorithms	Hyperparameters	Default value	Range searched	Selected value for Hs	Selected value for Tp	Selected value for Dir
	Neural networks Gradient Boosting Number of trees No. of units in hidden layer Activation function sigmoid 13 (2 × h + 1) {1-40} {relu,sigmoid,tanh} Learning rate 0.001 {0.0001-0.1} Epochs 30 {10,30,50,100,150} Batch size 32 {16,32,64,128} 100 {100-2000} trees Learning rate 0.3 {0.0001-0.3} Max depth 6 {1-20} Minimal node size 1 {1-15} Subsample 1 {0.5-1} Col sample 1 {0.5-1} Random forest Number of trees 500 {100,200,500,800,1000} Mtry 2 ( √ h) {2-6}	26 sigmoid 0.021 50 32 560 0.072 14 7 0.57 0.99 1000 2	20 sigmoid 0.016 100 64 1150 0.028 20 1 0.82 0.85 1000 2	40 relu 0.005 150 64 1990 0.069 20 1 0.79 0.9 1000 2

Table 3 :

 3 Statistical metrics for the three variables of interest after the hyperparameter tuning."Ann" stands for artificial neural networks, "Rf" for random forest and "Gb" for gradient boosting tree.

	Hs					Tp				θ p
	Numerical	Ann	Rf	Gb	Numerical	Ann	Rf	Gb	Numerical	Ann	Rf	Gb
	model	Corr.	Corr.	Corr.	model	Corr.	Corr.	Corr.	model	Corr.	Corr.	Corr.
						Computed with all data
	Biais	-0.201 0.026 -0.002 -0.001	0.712 0.007 -0.022 0.003	-0.249 0.790 0.979 0.714
	RMSE	0.399 0.300 0.246 0.240	1.839 1.553 1.269 1.231	13.803 12.07 9.646 9.501
	SI	0.166 0.144 0.118 0.116	0.153 0.140 0.114 0.111	0.045	0.04 0.032 0.031
	Cor	0.954 0.964 0.976 0.977	0.78 .817 0.882 0.889	0.330	0.36 0.461 0.421
	Biais	-0.536 -0.117 -0.120 -0.099	Computed with data where H s > 3m 0.683 0.032 -0.084 -0.051	1.925 -1.114 0.062 0.056
	RMSE	0.766 0.495 0.417 0.404	1.348 1.064 0.950 0.943	7.351 5.820 4.412 4.365
	SI	0.133 0.117 0.097 0.095	0.089 0.081 0.072 0.072	0.024 0.019 0.015 0.015
	Cor	0.818 0.861 0.903 0.908	0.832 0.856 0.888 0.889	0.589 0.609 0.793 0.793
	lower after the tuning.

Table 4 :

 4 Statistical metrics of the R2% calculated with Stockdon's formula. These results are obtained by taking the R2% computed with buoy data as reference. "Ann" stands for artificial neural networks, "Rf" for random forest and "Gb" for gradient boosting tree.

		Numerical	Ann	Rf	Gb
		model	Corrected	Corrected	Corrected
		Computed with all data
	Biais	0.003	0.019	0.002	0.003
	RMSE	0.223	0.209	0.175	0.172
	SI	0.145	0.136	0.114	0.112
	Cor	0.943	0.950	0.965	0.966
		Computed with data where Hs > 3m
	Biais	-0.042	-0.030	-0.054	-0.042
	RMSE	0.313	0.282	0.248	0.242
	SI	0.119	0.108	0.093	0.092
	Cor	0.854	0.862	0.897	0.901

tree algorithm decreases the RMSE of the extreme wave run-up by 22% (for all data and data where H s > 3m). Random forest shows almost the same reduction of RMSE values: 21.5% for all data and 20.7% for data where H s > 3m. The correction obtained by neural networks is less efficient: it reduces the RMSE computed with all data and data where H s > 3m by 6.2 and 9.9% respectively.

Table 1 .

 1 CNN performances for both sites. Best models are in bold font.

			(a) Biarritz
	CNN	Training Time (min)	Epochs	Time per Epoch (s)	Precision Recall	F 1 -Score
	Baseline Custom CNN AlexNet VGG16 Inception v3 Class weights Custom CNN AlexNet VGG16 Inception v3 Oversampling Custom CNN AlexNet VGG16 VGG16 Transfer Inception v3 Inception v3 Transfer	16.1 17.2 81.4 40.0 4.6 3.8 43.4 23.2 10.6 11.8 69.6 49.9 59.6 34.5	100 100 89 69 28 21 46 39 26 27 28 20 38 21	9.7 10.3 54.9 34.8 9.9 10.9 56.6 35.8 24.6 26.1 149.1 149.6 94.1 98.6	/ / / 0.721 / 0.568 0.574 0.563 0.642 0.716 0.783 0.869 0.679 0.777	0.333 0.333 0.481 0.714 0.603 0.777 0.832 0.798 0.880 0.885 0.851 0.865 0.767 0.786	0.328 0.328 0.476 0.713 0.474 0.609 0.645 0.631 0.718 0.777 0.813 0.866 0.717 0.780
			(b) Zarautz
	CNN	Training Time (min)	Epochs	Time per Epoch (s)	Precision Recall F 1 -Score
	Baseline Custom CNN AlexNet VGG16 Inception v3 Class weights Custom CNN AlexNet VGG16 Inception v3 Oversampling Custom CNN AlexNet VGG16 VGG16 Transfer Inception v3 Inception v3 Transfer	22.5 24.0 202.1 108.7 11.6 22.7 81.9 89.3 38.7 22.6 146.6 86.5 97.7 65.3	49 48 72 64 26 45 30 53 36 19 22 13 24 16	27.5 30.0 168.4 101.9 26.7 30.3 163.7 101.1 64.5 71.3 399.8 399.1 244.2 245.0	/ / / / 0.666 0.671 0.680 0.654 0.769 0.756 0.775 0.897 0.777 0.869	0.637 0.628 0.635 0.630 0.846 0.817 0.844 0.838 0.804 0.797 0.812 0.834 0.801 0.835	0.616 0.616 0.617 0.614 0.720 0.716 0.732 0.710 0.783 0.775 0.792 0.858 0.784 0.849
	4.1.2. Class Imbalance

Table 2 .

 2 Confusion matrices obtained by the best models for both sites.

(a)

Biarritz (best model: OV VGG16 Transfer)

			Swash	Predicted Collision	Overwash
	Observed	Swash Collision Overwash	1576 4 1	7 34 2	0 2 9
		(b)

Zarautz (best model: OV VGG16 Transfer) Predicted Swash Collision Overwash

	Observed	Swash Collision Overwash	4265 13 0	40 617 25	0 8 30

Table 3 .

 3 Classification results after correcting the misclassification in the test and validation sets.

	CNN	Time (

min) Epochs Time per Epoch (s) Precision Recall F1 Score

	Biarritz Best model before corr. Best model after corr. Zarautz Best model before corr. Best model after corr.	49.9 60 86.5 88.5	20 24 13 13	149.6 150 399.1 408.5	0.869 0.895 0.897 0.917	0.865 0.833 0.834 0.859	0.866 0.860 0.858 0.883

Table 4 .

 4 Performances of CNN learning from scratch, pre-trained with ImageNet, or pre-trained with the other site for Biarritz and Zarautz. "OV" stands for oversampling.

	CNN	Time (min) Epochs	Time per Epoch (s)	Precision Recall F 1 -Score
	Biarritz VGG16 (OV) VGG16 (OV) Pretraining with ImageNet VGG16 (OV) Pretraining with Zarautz data Zarautz VGG16 (OV) VGG16 (OV) Pretraining with ImageNet VGG16 (OV) Pretraining with Biarritz data	69.6 49.9 47 81.9 86.5 92	28 20 19 30 13 14	149.1 149.6 148.4 163.7 399.1 394.2	0.783 0.869 0.826 0.680 0.897 0.909	0.851 0.865 0.832 0.844 0.834 0.867	0.813 0.866 0.823 0.732 0.858 0.885

Table 1 :

 1 Table summarizing all the data available for the site of Biarritz.

	Data source	Variables	Period covered	Frequency
	Wave buoy Reanalysis and forecast wave models improved Weather station	Measured wave characteristics: H s , T p , θ p Simulated wave characteristics: H s , T p , θ p Meteorological conditions: Atm. Pressure, Wind speed and direction	11/2009-04/2020 30 minutes 01/1993-05/2021 Hourly 01/1993-11/2020 Hourly
	Tide gauge	Astronomical tide predicted with harmonical analysis.	01/1993-05/2021 Hourly
	Tide gauge	Atm. Surge extracted from measurements	05/2011-06/2020 Hourly
	Video monitoring station Storm impact regime extracted from images 03/2017-03/2021 15 minutes

Table 2 :

 2 Out of sample performances and optimal hyperparameters for the different models for the prediction of atmospheric surge. Results from 5-fold cross validation.

		RMSE SE RMSE
	Random forest	0.0565	0.0004
	Gradient boosting trees	0.0605	0.0004
	Shallow neural networks 0.1094	0.0007

Table 3 :

 3 Out of sample performances and optimal hyperparameters for the different models for the storm impact model. Results from 5-fold cross validation

		Kappa SE Kappa
	Random forest	0.6788 0.0265
	Gradient boosting trees 0.5648 0.0226
	Multinomial model	0.4734 0.0153

Table 4 :

 4 Performance of the best model (Random forest) on the training set

		Swash Impact Overwash
	Swash	8143	19	1
	Impact	69	131	3
	Overwash	4	11	39

Table 5 :

 5 Descriptive and predictive performances of the two BNs.

	Descriptive performance (training set)
		Precision Recall F1-score
	BN based on obs. data	0.743 0.709	0.72
	BN based on obs. and pred data	0.871	0.88	0.875
	Predictive performance (test set)
		Precision Recall F1-score
	BN based on obs. data	0.594 0.673	0.628
	BN based on obs. and pred data	0.63 0.789	0.691

Table 6 :

 6 Confusion matrices on the test set (predictive accuracy) obtained by both Bayesian networks.

	(a) BN based on observational data
			Predicted Swash Impact Overwash
	Observed	Swash Impact Overwash	1209 10 3	24 15 1	1 6 5
	(b) BN based on observational and predicted data
			Predicted Swash Impact Overwash
	Observed	Swash Impact Overwash	1205 10 0	23 19 2	6 2 7

Table A .

 A 2: Amplitude and phase of the first 10 harmonics estimated with the

	TideHarmonics package
	Harmonic components Amplitude Phase
	M2	1.337	291.5
	S2	0.468	306.1
	N2	0.283	281.9
	K2	0.132	301.6
	O1	0.071	69.7
	K1	0.063	157.4
	nu2	0.053	281.9
	mu2	0.050	268.8
	2N2	0.041	270.6
	Sa	0.040	208.9
	used to predict the astronomical tide (without the meteorological effect). The
	root mean square error (RMSE) computed on the test data is 10.4 centimetres.
	This result is consistent with the literature, for example Makarynskyy et al.
	(2004) obtained RMSE values between 8 and 15 centimetres for different neural
	networks architectures. An example of one month forecast with the harmonic
	analysis method is showed in the Figure A.2.

Table 1

 1 Mean bias and associated standard errors of the estimated parameters obtained by several regression methods for the exponential variance function. Results based on 500 replications.σ i = e γ|x Ti β| , γ = 0.02 , AR(1), α = 0.5

		λ = 0%		λ = 5%	λ = 10%
		β0	β1	β0	β1	β0	β1
	Least square	0.02(0.01) -0.02(0.01)	-0.72(0.01) 2.42(0.01)	-1.32(0.01) 4.62(0.02)
	M-estimation	0.02(0.01) -0.02(0.01)	-0.04(0.01) 0.12(0.01)	-0.22(0.01) 0.74(0.01)
	Generalized least square 0.02(0.01) -0.02(0.01)	-0.21(0.01) 1.09(0.01)	-0.32(0.01)	1.9(0.02)
	Weighted M-estimation	0.02(0.01) -0.02(0.01)	0.1(0.01) -0.38(0.01)	0.13(0.01) -0.51(0.01)
	Proposed method	0.02(0.01) -0.02(0.01)	-0.01(0.01) -0.18(0.01)	-0.08(0.01) -0.07(0.02)
	c	2.11		0.93		0.74
		σ i = e γ|x T i β| , γ = 0.02 , AR(2), α 1 = 0.5, α 2 = 0.2
		λ = 0%		λ = 5%	λ = 10%
		β0	β1	β0	β1	β0	β1
	Least square	-0.01(0.01) -0.02(0.01)	-0.69(0.01) 2.37(0.01)	-1.26(0.01) 4.53(0.02)
	M-estimation	-0.01(0.01) -0.02(0.01)	-0.04(0.01) 0.12(0.01)	-0.23(0.01) 0.76(0.02)
	Generalized least square -0.01(0.01) -0.02(0.01)	-0.22(0.01) 1.18(0.02)	-0.36(0.01) 2.11(0.02)
	Weighted M-estimation -0.01(0.01) -0.02(0.01)	0.1(0.01)	-0.4(0.01)	0.12(0.01) -0.51(0.02)
	Proposed method	-0.01(0.01) -0.02(0.01)	-0.06(0.01) -0.17(0.01)	-0.15(0.01) -0.07(0.02)
	c	2.11		0.91		0.72

Table 2

 2 Parameter estimates (β), their standard errors and z-values for the least squares and proposed method. The critical region of the significance test is of the full model |z| > 1.96.

		Least squares method	Proposed method with ĉ = 1.5 and σi = 1.03e 0.131|x i β|
		Estimate Std.Error Z-value	Estimate Std.Error	Z-value
	Intercept	9.378	0.157	59.765	8.784	0.018	493.763
	Dew point	3.214	0.122	26.446	2.953	0.011	268.731
	Temperature	-3.978	0.148	-26.845	-3.416	0.013	-260.880
	Pressure	-0.882	0.149	-5.909	-0.877	0.011	-77.211
	Cbwd NE	-1.455	0.104	-13.952	-1.117	0.021	-52.428
	Cbwd NW	-1.716	0.099	-17.284	-1.400	0.018	-77.689
	Cbwd SE	0.512	0.081	6.306	0.666	0.019	35.685
	Iws	-0.577	0.058	-9.888	-0.410	0.005	-89.071
	Is	-0.086	0.054	-1.575	-0.027	0.008	-3.475
	Ir	-0.453	0.039	-11.677	-0.380	0.004	-88.028
	Policy	0.362	0.343	1.058	0.222	0.023	9.464
	Time Policy Cos 3 years cycle	-7 × 10 -4 5 × 10 -4 -0.340 0.157	-1.373 -2.164	-6.10 -4 3 × 10 -5 -0.238 0.010	-18.052 -22.878
	Sin 3 years cycle	0.547	0.174	3.155	0.509	0.010	48.993
	Cos 2 years cycle	0.206	0.133	1.549	0.345	0.008	43.472
	Sin 2 years cycle	0.084	0.123	0.683	0.104	0.009	11.361
	Cos 6 months cycle	-0.342	0.123	-2.776	-0.355	0.008	-42.542
	Sin 6 months cycle	-0.365	0.119	-3.069	-0.425	0.008	-51.631
	Cos 4 months cycle	-0.184	0.120	-1.539	-0.038	0.008	-4.722
	Sin 4 months cycle	-0.035	0.121	-0.287	-0.084	0.008	-10.509
	Cos 3 months cycle	0.032	0.120	0.263	0.047	0.008	5.947
	Sin 3 months cycle	0.025	0.119	0.212	-0.005	0.008	-0.651
	lag1	/	/	/	0.746	0.005	162.733
	lag2	/	/	/	0.102	0.004	23.260

Interactive web application presenting the work: https://aureliencallens. shinyapps.io/application/

website https://rda.ucar.edu

Remote Sens. 2021,13, 1933. https://doi.org/10.3390/rs13101933 https://www.mdpi.com/journal/remotesensing

Acknowledgements

Acknowledgments

Funding was provided by the Energy Environment Solutions (E2S-UPPA) consortium and the BIGCEES project from E2S-UPPA ("Big model and Big data in Computational Ecology and Environmental Sciences"). The authors would like to thank the French national meteorological service "MeteoFrance" and Copernicus Marine Environment Monitoring Service for providing data.

Reproducibility

Meteorological data used in this article are private and can not be provided by the authors. However, the R code to perform the analysis and an example of data assimilation on wave forecast data (used in operational) are provided in this Github repository. Acknowledgments: The authors gratefully acknowledge European POCTEFA Program funding under the research project MARLIT EFA344/19 and the complementary program of OCA. The authors would like to thank Zarautz Town Council and the Directorate of Emergencies and Meteorology of the Basque Government for their continuous support of Zarautz's video station.

Acknowledgments

Funding was provided by the Energy Environment Solutions (E2S-UPPA) consortium and the BIGCEES project from E2S-UPPA ("Big model and Big data in Computational Ecology and Environmental Sciences"). The authors would like to thank the French national meteorological service "MeteoFrance" and Copernicus Marine Environment Monitoring Service for providing data.

Acknowledgements

This research was partially funded by the Australian Research Council project (DP160104292).

Funding: Funding was provided by the Energy Environment Solutions (E2S UPPA) consortium and the BIGCEES project from E2S-UPPA ("Big model and Big data in Computational Ecology and Environmental Sciences").

Data Availability Statement:

The images used to train the CNN are not publicly available due to the large size of the files. SIAME laboratory owns the Biarritz images data; they can be provided upon request. Images of Zarautz's site used in this work are from Azti and the Zarautz Town Council (shared). They can also be provided upon request. The python scripts and weights of the best CNN are available here: https://github.com/AurelienCallens/CNN_Timestacks, accessed on 1 May 2021. The web application used to label the images of Zarautz site (R programming language) is available here: https://github.com/AurelienCallens/Shiny_Classifier, accessed on 1 May 2021.

Errors of the numerical model for θ p depending on significant wave height measured at the buoy.

Abbreviations

The following abbreviations are used in this manuscript: The inception model was imported with Keras with the following function: keras.applications.InceptionV3().

The architecture is not displayed due to readability; the reader is referred to the original work of Szegedy et al. [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF] for more details. 2020), we tested random forest and gradient boosting tree to improve the numerical wave model predictions. With 5-fold cross validation on the data, we perform hyperparameter tuning and we choose the best generalizing model for each parameter.

Appendix A. 1. Reanalysis 19931. Reanalysis -2020 For this data assimilation task, random forest was the best generalizing algorithm for each of the 3 wave parameters. The metrics with and without the error prediction method are given in the table below. They are computed on all data and on data where there are stormy conditions (H s above 3 meters) because the wave model is known to greatly underestimate wave parameters in these conditions. For this data assimilation task, random forest was the best generalizing algorithm for each of the 3 wave parameters. The metrics with and without the error prediction method are given in the table below. They are computed on all data and on data where there are stormy conditions (H s above 3 meters) because the wave model is known to greatly underestimate wave parameters in these conditions.

Optimal hyperparameters

Random forest Ntrees = 500, mtry = 5 Gradient boosting trees Ntrees = 446, mtry = 6, lr = 0.011 Shallow neural networks Hidden units = 19, dropout = 0.232

Ntrees represents the number of trees, mtry represents the number of feature selected when forming each split in a single tree, lr designates the learning rate, hidden units corresponds to the number of neurons in the hidden layers and dropout represents the drop out rate of the neurons in the hidden layer.

Appendix B.2. Storm impact model

Table A4: Optimal hyperparameters found with the 5-fold cross validation.

Optimal hyperparameters

Random forest

Ntrees = 263, mtry = 3 Gradient boosting trees Ntrees = 428, mtry = 3, lr = 0.038 Multinomial model Penalty= 2.05e-10

Ntrees represents the number of trees, mtry represents the number of feature selected when forming each split in a single tree, lr designates the learning rate, Penalty corresponds to the penalty of the multinomial model fitted by the R package glmnet.

Appendix C. Training the bayesian networks

Appendix C.1. Discretization of the variables

Table A5: Discretization thresholds chosen for the explanatory variables.

Variables

Number of class Ranges for class Metric

Significant wave height 5 [0,1.5]; [1.5,3]; [3,4.5]; [4.5,6]; [6,10] m

Wave period at peak 5 [0,7]; [7,9]; [9,12]; [12,15]; [15,22] s

Wave direction 3 [0,300]; [300,305]; [305,360] degree Tide 6 [0,2]; [2,3]; [3,3.5]; [3.5,4]; [4,4.5]; [4.5,6 in this thesis (chapter 3 and 5) to obtain satisfactory results. However, we saw the limitations of these methods in chapter 5, where we have only a few observations of storm events. It would be interesting to extend the work of this chapter by combining the observational database with a synthetic database created in the same spirit as proposed by Poelhekke et al. (2016) which we have exploited in Morichon et al. (2018) for the Grande Plage Biarritz. In the study of Poelhekke et al. (2016), they employed a statistical method (copula), which capture the dependencies between the storm characteristic variables, to simulate synthetic storm events. The impact of these storm events were then simulated by a process-based model. This methodology combining a statistical method and a process-based model could be the solution to improve the storm impact model presented in last chapter. In any case, this limitation highlights the need to maintain existing monitoring networks and to install new ones in order to train better statistical learning models for the prediction of coastal risks.

Temporal aspect

In this research work, the temporal aspect has not been treated. It could be interesting to integrate this temporal aspect particularly in the improvement of wave forecast made by spectral wave model. Indeed, it is reasonable to think that the errors of the spectral wave model are temporally correlated.

To include this temporal aspect, we could use errors of the previous time step as explanatory variables and we could also employ SLM that are adapted for time series modeling such as recurrent neural networks [START_REF] Zhang | Ocean wave height series prediction with numerical long short-term memory[END_REF] or random forest adapted to time series [START_REF] Goehry | Random forests for time series[END_REF]. The temporal aspect could also be included in the storm impact model by employing dynamic bayesian networks which are bayesian networks including the concept of time [START_REF] Murphy | Dynamic bayesian networks: Representation, inference and learning[END_REF].

maximum number implemented in the R package TideHarmonics. Tide harmonic method was chosen over neural networks for the modelling of the tide level because once the model is fitted, past observations are not needed to make a prediction unlike neural networks. This is convenient as the data from Socoa tide gauge contains a lot of missing data.

Results and discussion

The estimated amplitudes and phases for the 10 first harmonic components of Socoa are presented in the table A.2. The 2 harmonic components with the highest amplitudes are M2 and S2 which are semi-diurnal components. This indicates that the type of the tide in the area of Socoa is semi-diurnal [START_REF] Nicolle | Modélisation Des Marées et Des Surcotes Dans Les Pertuis Charentais[END_REF].

Once all the harmonic amplitudes and phases are estimated, they can be Noname manuscript No.

(will be inserted by the editor)

Robust estimation procedure for autoregressive models with heterogeneity

Robust estimation procedure for autoregressive models with heterogeneity 8 squares method (gls function from the nlme package),M-estimation with c = 1.345 (rlm function from the MASS package), the weighted M-estimation (whm function from rlmDataDriven package) with c = 1.345 and our datadriven method (rlmDD het function from rlmDataDriven package).

For one simulation, we generate a multivariate normal dataset (n = 500) using the model (1). In our case, x T i β = β 0 + β 1 x i1 where we fix the value of β 0 and β 1 to 10 and x i1 comes from a uniform distribution on (0, 1). For σ i , we test two functions: the power function σ i = |x T i β| γ with γ = 0.2 and the exponential function σ i = e γ|x T i β| with γ = 0.01. For the term i , we consider two cases: (i) an autoregressive process of order 1 with α = 0.5 and (ii) an autoregressive process of order 2 with α 1 = 0.5, α 2 = 0.2. These processes can be written as follows: i = p j=1 α j i-j + ξ i where ξ i are independent and normally distributed errors following a standard normal distribution N (0, 1). In order to add outliers, the data is randomly contaminated by adding a value sampled from a uniform distribution on (0, 10). Several contamination rates are considered: λ = 0%, 5%, 10%.

Table 1 shows the means and the standard errors of the estimate bias computed accross all the simulations for the exponential function. The results obtained with the power function are included in the online supplementary material (Table S-1).

For both variance functions, average value of the data-dependent tuning constant decreases as contamination rate increases. This tuning constant has the expected behaviour: as the proportion of outliers becomes larger, more values should be considered as outliers therefore the tuning constant is smaller.

The bias and associated standard errors both increase with the proportion of outliers. However it is less important for robust methods, especially for the proposed method. Indeed, when data are contaminated, the method yields estimators with lower mean bias and associated standard errors.

This method also provides estimates with lowest MSE in all the cases when contamination is present. Differences in efficiency with the other methods augment with the contamination rate, reaching MSE values fives fold lower than the ones obtained by non robust methods (Figure 1). The same figure for the power variance function is available in the supplementary material (Figure S-1).