Atmospheric storm surge modeling methodology along the French (Atlantic and English Channel) coast

Storm surge modeling and forecast are the key issues for coastal risk early warning systems. As a general objective, this study aims at improving high-frequency storm surge variations modeling within the PREVIMER system (www.previmer.org), along the French Atlantic and English Channel coasts. The paper focuses on (1) sea surface drag parameterization and (2) uncertainties induced by the meteorological data quality. The modeling is based on the shallow-water version of the model for applications at regional scale (MARS), with a 2-km spatial resolution. The model computes together tide and surge, allowing properly taking into account tide-surge interactions. To select the most appropriate parameterization for the study area, a sensitivity analysis on sea surface drag parameterizations is done, based on comparisons of modeled storm surges (extracted with a tidal component analysis) with four tidal gauges, during four storm events, and over about 7.5 years, where the observed water level is processed in the same way as the modeling results. The tested drag parameterizations are a constant one, as reported by Moon et al. (J Atmos Sci 61: 2321–2333, 2007), Makin (Bound-Layer Meteorol 115: 169–176, 2005), and Charnock (J Roy Meteor Soc 81: 639–640, 1955). Charnock’s parameterization, either constant with high value (0.022) or relying on a full statistical description of the sea state, enables to improve storm surges forecast with peak errors 10 cm smaller than those computed with the other drag coefficient formulations. The impact of the meteorological forcing quality is evaluated over January 2012 from the comparison between surges modeled with different meteorological data (ARPEGE, ARPEGE High Resolution and AROME) and observations. For event time scale, storm surge computation is highly improved with ARPEGE High Resolution data. For month time scale, statistics of model accuracy are less sensitive to the choice of meteorological forcing. As a conclusion, the Charnock’s parameterization is advised to model storm surges on the French Atlantic and English Channel coasts, whereas the quality requirements regarding meteorological inputs depend on the time scale of interest. Within the PREVIMER system, aiming at forecasting events, ARPEGE High Resolution data are used.


Introduction
The recent coastal flooding events like during hurricane Katrina of 2005 or the mid-latitude extra-tropical storm Xynthia that hit the central part of the Bay of Biscay in 2010, illustrate the need to provide realistic sea level forecast. The mean water level along the coast results from tide, surge, and wave setup. Most of the forecast systems provide water level including tide and surge components, based on different strategies taking tide-surge interaction or not into account (Idier et al. 2012a). Moreover, operational surge models whose spatial resolution is about a few kilometers use either a constant sea surface drag coefficient or those issued from Wu (1982), Smith and Banke (1975), or Charnock (1955) formulations without taking into account the sea state. For instance, for WAQUA-IN-Simona/DCSM98 (Verlaan et al. 2005), DaMSA (Büchmann et al. 2011), andPrevimar (Le Breton andPerherin 2009) systems, root mean square errors are about 5 to 15 cm for period covering a few days or several years. In France, in addition to Previmar, a forecast system has been developed (PREVIMER www.previmer.org) (Lecornu and De Roeck 2009) providing water levels and surges along the French Atlantic and English Channel coast. The surge modeling strategy relies on embedded models, with a total of seven models (Fig. 1). The rank 0 model covers North East Atlantic with a spatial resolution of 2,000 m and provides total water levels to the rank 1 model which covers the Bay of Biscay and the English Channel with a resolution of 700 m. The five rank 2 have a resolution of 250 m, whereas their boundary conditions is the sum of the storm surge computed from the rank 1 and tidal conditions (Le Roy and Simon 2003). This system provides water levels and storm surges forecast of 15-min temporal resolution. The first version of this forecast system (2008) was based on a model having a constant sea surface drag coefficient and using meteorological data from ARPEGE model of Meteo-France (Courtier et al. 1991(Courtier et al. , 1994. Recent studies like the modeling work of Bertin et al. (2012a) on the Xynthia storm surge highlight the necessity of including drag formulations better than a constant one in storm surge models. In addition, recent progress has been done on meteorological models including now orographic phenomena, as well as a better physical description of mesoscale processes (e.g., the AROME model, Seity et al. 2011), and providing higher resolution meteorological data.
The aim of the study is to improve the large scale model of the system (rank 0) in terms of storm surge reproduction skills, using the progresses done on drag coefficient parameterization, as well as meteorological models. The method is based on a sensitivity analysis of the model used in the PREVIMER system, regarding sea surface drag coefficient formulation and meteorological data. The approach focuses not only on storm events but also on multi-annual time span especially for the drag sensitivity analysis.
The paper is organized as follows. First, the model is presented (Sect. 2). Then, a sensitivity study to sea surface drag formulation is carried out leading to propose an improved storm surge model for the forecast system (Sect. 3). In Sect. 4, sensitivity to the meteorological data is tackled, using the improved model of Sect. 3. The results are discussed in Sect. 5, before the conclusion is drawn.

Description of the reference configuration
Water level computation is issued from high frequency sea level variations modeling with the shallow-water version of the model for applications at regional scale called MARS 2DH, developed by IFREMER (Lazure and Dumas 2008). The modeling approach followed in this study is structured around a reference configuration whose specificities are described below. The configuration under study is a 2-kmresolution one and covers the North East Atlantic (Fig. 2).
The bathymetry comes from the 1′ North West European shelf Operational Oceanographic System (NOOS, http:// www.noos.cc) bathymetry and 100-m resolution digital terrain models from IFREMER and SHOM (French Hydrographic and Oceanographic Service). The bottom friction τ ! b is parameterized with a quadratic formulation defined as where ρ is the density of the sea water, g is the acceleration due to gravity, u ! is the current, h is the depth, and k is the Strickler coefficient set equal to 35m 1/3 s −1 . The sea surface friction is also computed with a quadratic parameterization.
where ρ a is the air density (set to 1.25), C D is the constant drag coefficient equal to 0.0016, and U 10 ! is the 10-m wind.
As far as meteorological forcing are concerned, the reference configuration uses the 0.5°ARPEGE meteorological model (pressure and wind) (provided by Meteo-France, Courtier et al. 1991Courtier et al. , 1994, with a 6-h time step. Before 2009, it consists of four analyses a day, and since 2009, it consists of analysis at 0000 h and forecasts at 0600, 1200, and 1800 h. These meteorological fields are linearly interpolated over MARS model grid.
3 Sensitivity to sea surface drag parameterization

Methodology of comparison with observations
To identify the sea surface drag formulation that best reproduces storm surge along the studied area, a comparison of model outputs issued from the sensitivity study with tidal gauge data has been carried out over a period covering 7.5 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). The water level data are issued from the REFMAR database (http://refmar.shom.fr/en) which contains water level observations recorded by tidal gauges along the French coast. They consist in hourly data which result from the smoothing of 10-min data. Here, we focus on storm surges. Water level observations, after some treatment, can provide the so-called practical surge, obtained by subtracting the predicted tide from the observed water level. One technique to estimate surges is to make double runs (one with tide and meteorological forcing, one with tide only) and then obtain surge from their difference. However, the predicted tide is often provided by multi-decadal water level signal processing, whereas the modeling is done only on few events or years. Furthermore, tidal component analysis also includes radiational components which are related to the meteorology. Thus, in order to properly compare our model results to observations, the same processing is applied to the model and data, and over the same time span, the water level signal (issued either from the model or from the observation) is first processed by a tidal component analysis, providing a predicted tide, and then by subtraction from the total water level, the modeled and observed practical surges are obtained. The tidal analysis has been computed, thanks to the French Hydrographic Office software MAS (Simon 2007), which computes until 143 harmonic constituents and takes into account nodal corrections. The selected time span is the 2003-2010 period.
The storm surges computed from the different sea surface drag parameterizations are evaluated in relation to observations over four short representative storm events and over the 2003-2010 time span. It enables to estimate if the mean behavior of the modeled surge signal and the storm surge dynamics during strong meteorological events is well reproduced. The four specific storm events (Fig. 3)   The sensitivity analysis focuses on four locations of interest ( Fig. 2), Dunkerque, Saint-Malo, Le Conquet, and La Rochelle, and is based on error indicators (listed in Table 1) to conclude about the most suitable sea surface drag formulation for storm surge modeling.
Among the error indicators, MEAN, root mean square error (RMSE), COST, BIAS, and EFF give information about the similarity between the mean behavior of the surge signals issued from modeling and observations, whereas max error (T=½year, T=1 year, T=2 years), peak error, and phase shift reveal the ability of the model to reproduce high storm surges. Error (T=½year, T=1 year, T=2 years) is based on the following extreme statistics analysis applied to the surge time series (modeled and observed) over the 7.5 years. From the peak over threshold (POT) method, independent events are spotted when the surge is higher than the first centile of the surge series. Events separated by at least 3 days are considered as independent. Then, a generalized Pareto distribution (GPD) law is adjusted on the surge sampling from the maximum likelihood method (Coles 2001); it gives the probability of surge level to be smaller than the return surge h (P(X≤h) with X as the independent surges. The return period T associated to the surge level h is defined as =1/(1−GPD(h))n, with n as the number of storm events per year. Then, the return surge is h= GPD −1 (1−1/nT).

Sea surface drag formulations and comparison
The sea surface drag coefficient formulation plays a key role in storm surge modeling (Nicolle et al. 2009;Sheng et al. 2010). According to shallow water equations which are the basis of the surge model and under some assumptions (1D uniform steady flow over a flat bed with wind directed to the coast, neglecting Coriolis effect), it is possible to analytically Only used for storm event analysis Root mean square error (RMSE) Nash and Sutcliffe (1970). EEF is close to 1 for model results close to observations Peak error Phase shift is the difference between times of the maximum storm surge in observations and model outputs over the studied period Errors are computed on the instantaneous practical storm surge. Where Di and M i represent respectively the observed and modeled storm surge, N is the number of data, and i stands for the instantaneous value show that the storm surge is increasing with the sea surface wind stress. Equation (2) shows that, for a constant wind, the wind stress increases with the drag coefficient. Thus, as a first approximation, for a constant wind, the surge should increase with the drag coefficient. This behavior has been shown in many studies, as for instance Moon et al. (2009) and Zweers et al. (2011). For a deeper analysis of the relation between wind, wind shear stress, drag coefficient, and skew storm surges, see the study of Zweers et al. (2011). The first level of drag coefficient is the constant value. Going more in details, the sea surface drag coefficient is a function of wind speed, wave state, and atmospheric stability (Monin and Obukhov 1954;Charnock 1955;Large and Pond 1981;Smith et al. 1992). However, the parameterizations employed for this coefficient are often only wind-dependent in ocean modeling. Recent results issued from in situ measurements, lab experiments, or theoretical works (Alamaro et al. 2002;Powell et al. 2003;Donelan et al. 2004;Moon et al. 2004a, b;Makin 2005) show that the sea surface drag coefficient is weaker under strong wind (∼30 m/s) conditions and could even be a decreasing function of the wind speed. It means that sea surface drag coefficient could be sizeably dependent to the wave state (Toba et al. 1990;Smith et al. 1992;Johnson et al. 1998;Drennan et al. 2003). Indeed, the aerodynamic roughness length associated with the surface wavefield (z 0 ) impacts atmospheric circulation and storm surges (Mastenbroek et al. 1993;Saetra et al. 2007;Nicolle et al. 2009). Assuming a logarithmic variation of wind speed with height (z) above the surface, sea surface drag coefficient is linked to the surface roughness, for 10-m wind, when z=10m through the following equation: where u * is the friction velocity and κ is the von Karman's constant.
A first class of drag formulation is the 10-m wind-dependent parameterizations (Moon et al. 2007;Makin 2005). Moon's empirical formulation enables to take air-sea momentum transfers for strong winds into account with more accuracy. The surface roughness is linked to 10-m wind using the law Therefore, using the Eq. (3), the sea surface drag coefficient can be expressed at 10 m, as a function of 10-m wind. Concerning Makin's formulation, the expression of the surface roughness is where w ¼ min 1; acr = κuÃ Þ; a cr ¼ 0:64m=s; c l ¼ 10 and c z0 ¼ 0:01 ð . The friction velocity and consequently the sea surface drag coefficient can be iteratively computed by replacing z 0 in Eq. (3) with this formulation.
A second class of drag formulation is based on Charnock's (1955) wave-and wind-dependent formulation. According to Charnock (1955), the surface roughness is expressed by where α c is the Charnock dimensionless parameter. Then, the drag coefficient can be computed from this equation and Eq. (3). In the present study, α c is constant or variable (α cvar ) and, in this case, issued from the IOWAGA modeling system (Rascle and Ardhuin 2013) which is based on the wave model WAVEWATCH III® (Tolman 2008;Ardhuin et al. 2010;Tolman et al. 2013). Within the wave model, α c comes from Eq. (6), thanks to the iterative computation of the friction velocity using the wave growth rate following the parameterization of Janssen et al. (1994) as used at European Centre for Medium-Range Weather Forecasts (ECMWF). An important difference, though, is that the parameterization for the wave breaking and generation terms are different from those of ECMWF (Bidlot 2012), in particular because Ardhuin et al. (2010) include a reduction of the wavesupported stress in the high frequencies to represent a sheltering of short waves by longer waves. Although this parameterization leads to better short wave properties, it removes most of the wave-induced variability in the Charnock coefficient (Rascle and Ardhuin 2013). We thus expect that our results are much less sensitive to the wave age than using other parameterizations (e.g., Mastenbroek et al. 1993). The IOWAGA system provides α cvar with a 0.5°resolution and a time sample of 3 h. The meteorological forcing used to run the wave model comes from Climate Forecast System Reanalysis (CFSR) (Saha et al. 2010) up to 2006 and then from the ECMWF up to 2010.
Using these formulations in our model over the year 2008, we obtain different values for the drag coefficient depending on wind and wave conditions (Fig. 4). Moon and Makin drag coefficients are plotted, along with Charnock drag coefficients over the year 2008, issued either from ARPEGE winds and α c =0.014 (Ardhuin, pers. communication) or computed from the IOWAGA modeling system (α cvar ). For wind ranges of 2 m/s, C d (α c =cst) and C d (α cvar ) over the year 2008 and the whole area of interest (at each ARPEGE grid point) are represented by a point (green or blue, respectively), whereas continuous lines (green and blue) represent their mean value. Figure 4 highlights the differences between Moon, Makin, and Charnock's formulations for wind speed over 10 m/s. Moreover, for wind speed ranging between 20 and 33 m/s, C d (α cvar ) gives drag coefficients 1.5 times larger than the Moon or Makin formulations. The behavior of C d (α c =0.01) (not shown here) is almost the same as those of Makin drag coefficient for wind speed under 30 m/s (Makin 2005). For wind speed between 30 and 33 m/s, drag coefficient values issued from Moon et al. (2007) and Makin (2005) are close. For stronger wind, curves follow opposite trends: whereas Moon et al. (2007) relation keeps increasing and the drag coefficient issued from Makin (2005) formulation decreases. The comparison between C d (α c =0.014) and C d (α c =0.022) shows that a higher α c constant values gives results closer to those with α cvar . Nevertheless, Fig. 4 also highlights the fact that setting α c to a high constant value (0.022) for drag coefficient computation is not enough to obtain the same values as C d (α cvar ). For wind speed under 21 m/s, C d (α c = 0.022) is larger than C d (α cvar ), whereas for wind speed above 21 m/s, it is the opposite trend. In the next two sections, we analyze the modeled surge sensitivity to the wind-dependent drag formulations and then to the wind-and wave-dependent formulations.

Pluri-annual scale analysis
The first step of the sensitivity study on sea surface drag parameterizations focuses on those which are wind-dependent: Moon et al. (2007), Makin (2005), and Charnock (1955) with a constant α c (equal to 0.014 or 0.022). The results issued from MARS simulations with these three formulations are compared to those issued from the reference simulation and observations. Error indicators (Table 1)  Rochelle and about 5 cm at the Le Conquet whatever the sea surface parameterization used. Similarly, for each tested parameterizations, the absolute value of bias is about 0.8 cm at Dunkerque, 1.5 cm at Saint-Malo, and smaller than 0.3 cm at Le Conquet and La Rochelle. The sensitivity of the surge signal to the different sea surface drag coefficient formulations is more obvious for error indicators relative to large storm surge (errors on the surges associated with ½-, 1-, and 2-year return periods, phase shift, and peak errors). The errors associated with ½-, 1-, and 2-year return period surges should be interpreted in light of the observations (see Table 2). On the whole, Charnock (with α c =0.022) and Moon formulations give the smallest errors on the surges associated with the different return periods, with errors ranging from 13 to 44 cm for observed surges associated with ½-, 1-, and 2-year return periods varying between 0.52 and 1.57 m. The same trend is noticeable on peak errors. Phase shifts are the highest at Le Conquet (∼3.4 h) and the smallest at Dunkerque (∼0.9 h). If we focus on each location, at Dunkerque (Fig. 5), Charnock (with α c =0.022) formulation reduces errors associated with ½-, 1-, and 2-year events up to 10 cm compared to Moon formulation. At Saint-Malo (Fig. 6), Moon and Charnock (with α c =0.014) formulations give similar results concerning errors on the surges associated with ½and 1-year return that are the most satisfactory. But for errors on the surges associated with 2-year events and peak error, Moon and Charnock (with α c =0.022) formulations prevail over the others. At Le Conquet (Fig. 7), Charnock (with α c = 0.022) formulation gives the smallest errors relative to high storm surge. It is followed by Moon and Charnock (with α c = The fact that the difference between the formulations is more obvious on the statistics relative to high surges than to those which describe the general trend of the surge signal is in agreement with the behavior of the drag coefficients as a function of wind. Indeed, for weak wind speed (about 10-15 m/s), drag coefficients are close to 0.0016, whereas for strong wind (>25 m/s), their deviations can be 1.5 times higher than the value 0.0016 (see Fig. 4). These differences between drag coefficients imply differences in terms of storm surges which are all the greater that the wind is strong. The long-term analysis shows that even if the results issued from the different sea surface drag formulations are similar, Charnock (with α c = 0.022) parameterization gives the more accurate ones. Errors on the storm surges associated with ½-, 1-, and 2-year return periods and peak errors are reduced with Charnock (with α c = 0.022) parameterization. The reason seems to be the higher values for drag coefficient issued from Charnock (with α c = 0.022) parameterization than Makin (see Fig. 4). As an example, the mean value of Charnock drag coefficient (with α c = 0.022) over a square area around Dunkerque with a 200-km extent and over 2003-2010 is 15 % higher than the mean value of Makin drag coefficient over the same geographical extent and period. Nevertheless, the example of Saint-Malo shows that for some location, Charnock (with α c =0.014) parameterization can give the more accurate results for high surge computation if we focus on errors associated with ½and 1-year events, which justify the use of a time and space variable α c .

Event scale analysis
The short-term analysis was performed for the same locations as the long-term one. It concerns the four storm events described in part 3.1. For each studied storm events, Charnock formulation with α c =0.022 gives the more satisfactory results compared to the observations (Figs. 9, 10, 11, and 12). As error indicators, we focus on maximal and peak errors which  (Fig. 13). The improvement of the storm surge modeling with Charnock formulation (with α c =0.022) is less obvious at Le Conquet and La Rochelle during Johanna storm (Fig. 10b, c) since the gaps between drag coefficients are weaker nearby these locations (Fig. 14). For Quentin storm event, the peak error is reduced by 22 cm at La Rochelle with C d (α c =0.022) regarding the reference configuration with C d =0.0016 (Fig. 11, Table 3). In terms of drag coefficient, Moon formulation gives a value 70 % higher nearby La Rochelle than the constant one (Fig. 15), which explains the improvement of the modeled storm surge peak's height. At Dunkerque and Saint-Malo, the contribution of the Charnock formulation (with α c =0.022) to storm surge modeling is not obvious. At Le Conquet, the decrease of the surge after February 10, 2009 is better reproduced by the model using a Charnock formulation (with α c =0.022) than the baseline (Fig. 11c). For Xynthia storm event (Fig. 12), Charnock formulation (with α c =0.022) reduces the peak error by 16 cm at La Rochelle. For the other locations, all the tested formulations give results very similar to the reference configuration. It is coherent with the locations where the gaps between Moon, Makin, and Charnock (with α c =0.014 or α c =0.022) coefficients and the value 0.0016 are the smallest (Fig. 16). This storm event analysis shows the added value of Charnock formulation with α c =0.022 in storm surge modeling, keeping in mind that it can induce some side effects by overestimating local drag coefficients. The sensitivity study to sea surface drag parameterization has revealed the relevance of the use of Charnock formulation (with α c =0.022) in storm surge modeling compared to Moon, Makin, and Charnock (with α c =0.014) formulations or to a constant drag coefficient, thanks to the long-term and short-term analysis. The error indicators computed for the short term analysis, which are not presented in the paper, showed that Charnock formulation with α c =0.022 reduces significantly only peak and maximal errors (up to 40 cm) at some locations (Dunkerque and La Rochelle) during special storm events (storm of November 9, 2007 and Quentin storm). Nevertheless, the RMS and mean errors only differ from a few centimeters whatever the formulation used. The statistical analysis realized on the high storm surges over 7.5 years has also proved the contribution of Charnock formulation (with α c =0.022) to storm surge modeling by reducing errors relative to high surges up to 19 cm at Dunkerque (error associated to 2-year return surge) compared to baseline. The added value of Charnock formulation (with α c =0.022) in surge modeling mostly concerns the surge peak reproduction.

Wind-and wave-dependent parameterization
To go one step further in this sensitivity study on sea surface drag coefficient, Charnock formulation with a variable α c has also been tested. This parameterization takes the sea surface roughness induced by waves into account. For instance, the analysis of the January 2008 period shows that the Charnock drag coefficient tends to increase with wave steepness, especially for large windspeed (not shown here).

Pluri-annual scale analysis
As for wind-dependent parameterizations, the error indicators relative to the mean storm surge signal, issued from multiannual comparisons between the modeled and the observed storm surge, reveal few differences between MARS simulations (Figs. 5, 6, 7, and 8). The improvement of the reproduction of the storm surge by Charnock formulation with a variable α c compared to the reference configuration (constant drag coefficient) and only wind-dependent parameterizations is also obvious as far as the error indicators relative to strong storm events are concerned. Actually, Charnock formulation with a variable α c gives on the whole, the smallest peak errors, and errors on the surges associated with the different return periods, compared to the reference configuration and also the configurations which use Moon or Makin parameterizations. Nevertheless, the comparison between simulations based on C d (α cvar ) and C d (α c =0.022) gives similar results. Differences between the error indicators associated to high surges concerning these simulations are smaller than 4 cm. Besides, at Dunkerque, using C d (α c =0.022) can even lead to a better match (improvement of 2 cm) between model and observations compared to the other formulations. However, even if results are bunched together (Figs. 5, 6, 7, and 8), the comparison between the simulation based on a Charnock formulation using a variable α c and those using a constant one has revealed the importance of taking into account concretely sea state in surge modeling. At Dunkerque, the errors on the surges associated with the different return periods are reduced respectively up to 19 and 15 cm, thanks to Charnock formulation with α c =0.022 and with a variable one (see Figs. 5, 6, 7, and 8) compared to baseline. At Saint-Malo, C d (α cvar ) reduces the errors relative to return surges up to 9 cm, whereas those with C d (α c =0.022) reduces them by 5 cm. At La Rochelle, the same trend is observed with a 6cm, resp. 3 cm, decrease respectively with the configuration with C d (α cvar ), resp. C d (α c =0.022), toward the reference configuration. Figure 4 illustrates these observations since the configuration with C d (α cvar ) gives the highest drag coefficients for wind speed over 21 m/s.

Event scale analysis
The short-term analysis does not confirm significatively the added value of a wind-and wave-dependent parameterization in storm surge modeling. For November 2007 and Johanna storms, as Charnock formulation with α c =0.022, the one with C d (α cvar ) enables to reduce maximal and peak errors by 40 cm at Dunkerque and by more than 20 cm at Saint-Malo (Figs. 9 and 10b and Table 3) compared to the reference configuration. The maps of the maximum difference between C d (α c =0.022), as well as C d (α cvar ), and the value 0.0016 show that the deviations can locally reach 80 % (Figs. 13 and 14). For Quentin storm, the peak error is reduced respectively by 26 cm, resp. 22 cm by using a configuration with C d (α cvar ), resp. C d (α c =0.022), at La Rochelle compared to reference (Fig. 11d). In terms of drag coefficient, C d (α cvar ) gives a maximal value about 0.003 nearby La Rochelle that is two times higher than the constant value of 0.0016 used in the  Table 3 and Figs. 9, 10, 11, and 12). From the error indicators related to high surge issued from the above long-term and short-term analysis and also the sensitivity study to wind-dependent formulations, it is concluded that Charnock formulation with a variable α c gives the more satisfactory results for the present storm surge modeling. However, the improvement in surge modeling is more due to a global increase of the Charnock coefficient than to its spatiotemporal variations, compared to Moon and Makin formulations or a constant drag coefficient. The similitudes between the maps of the maximum difference in percent between sea surface drag coefficient C d (α c =0.022) or C d (α cvar ) and the reference value (0.0016) over the studied storm events enable to anticipate this conclusion (see Figs. 13,14,15,and 16). Despite that, there are still some differences between the patterns described by these maps, which can lead to some differences in the computed surges outside the studied locations. Besides using a variable α c has more physical sense than a constant one; therefore, even if its added value is not so significant, it makes sense to integrate it in surge modeling, in order to avoid overestimation in local areas where the sea roughness is smaller, like for instance in protected environments. For these reasons, the drag formulation with wind-and wave-dependent Charnock coefficient is chosen for the storm surge forecast model within the PREVIMER system.

Meteorological data description
In order to investigate the sensitivity of the forecast model to time and space resolution of meteorological forcing, three different meteorological datasets from Meteo-France were  Table 4. Their extensions are presented in Fig. 17. Daily values at 0 h correspond to analysis; other values correspond to forecasts.
In order to estimate the meteorological forcing quality, meteorological model outputs (wind and pressure) were compared not only to local meteorological stations (Dunkerque, Saint-Malo, Le Conquet, and La Rochelle) but also to QuikSCAT scatterometer wind observations (Idier et al. 2012b). Here, only the comparison between meteorological data measured at Dunkerque station (hourly data provided by Meteo-France) and meteorological models is presented. Wind values correspond to mean winds over 10 min before the hour, measured at 10 m above the sea. Atmospheric pressure corresponds to the pressure reduced at sea level, measured every hour. These data have been compared with AROME, ARPEGE, and ARPEGE HR meteorological models in January 2012. Comparisons for atmospheric pressure, wind direction, and wind intensity from the 4th to the 8th of January 2013, during Andrea storm, are presented in Fig. 18. Atmospheric pressure and wind direction from the three models are very similar to measurements. Concerning wind intensity, AROME and ARPEGE are quite close to measurements but tend to underestimate wind intensity, whereas ARPEGE HR tends to overestimate wind intensity. In January 2012, for AROME, ARPEGE HR, and ARPEGE, the mean quadratic error is respectively 2, 2.9, and 1.8 m/s, and the bias is respectively −1.4, 2, and −0.5 m/s. The negative bias for ARPEGE HR confirms that the model tends to overestimate wind intensity of about 2 m/s.

Methodology
The model configuration under study is the optimized model presented in Sects. 2 and 3; its spatial resolution is 2 km (Fig. 1). The drag coefficient is wind-and wave-dependent; it is computed from the Charnock formulation. Charnock coefficients come from global WAVEWATCH III® model of IOWAGA project (Ardhuin et al. 2010) as described in Sect. 3.2.
Model experiments are done for the December 2011-January 2012 period. In the 15th and 16th of December 2011, Joachim storm crossed the France. As in Sect. 3, the analysis is done on the following four locations: Dunkerque, Saint-Malo, Le Conquet, and La Rochelle (Fig. 1).
The winds reached 32.2 m/s at Brest (near Le Conquet) and 30 m/s at La Rochelle. At the beginning of January 2012, two storms occurred: Illi storm, which crossed the north of France and the UK, on the 3rd of January and Andrea storm on the 5th of January.
The methodology to compute the modeled and measured surges is different from the methodology presented in Sect. 3, because in this case, a short operational run (2 months) is exploited instead of a long hindcast (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). Modeled surges are the differences between the simulation with both tide and meteorological forcing (wind and atmospheric pressure) and the simulation with only tide. This way, surges include also tide-surge interaction (Idier et al. 2012a). Measured surges have been computed from tide gauges. Water levels come from SHOM Sea Level Observation Network (RONIM, Poffa et al. 2012), and they have been analyzed with the same SHOM tidal software MAS (Simon 2007) as in Sect. 2, with the difference that analyzed time series start in 1996 for Dunkerque, 2003for Saint-Malo, 1992 for Le Conquet, and 1997 for La Rochelle. This software computes a maximum of 143 harmonic components (depending on the duration of the measurements: the longer is the measurement period, the more important is the number of harmonic constituents) and takes into account nodal corrections. Between 127 and 137, harmonic components were computed in each harbor. In order to have good quality harmonic components, only measurements from numeric tide gauge have been analyzed, which are generally periods long enough to ensure a good harmonic analysis (generally more than 10 years, often 30 years). From these precise harmonic components, predictions with MAS have been made in December 2011 and January 2012. Measured surges are the  -11, 2007 differences between water level measurements and tidal predictions. Despite all these precautions, there is often still a periodic signal in the measured surges. An example is presented in Fig. 19. This figure, extracted from PREVIMER website, shows surges in Saint-Malo that have been computed by French Hydrographic Office, which is responsible of the official tide prediction in France. Modeled surges are in blue (model with meteorological forcing minus model without meteorological forcing), and measured surges computed at tide gauge are in red. The oscillation range can exceed 70 cm. They are relatively permanent in Saint-Malo. The presence of this semidiurnal signal is sometimes due to the shift phase between measurement and tidal prediction, because of meteorological effects. The non-stationarity of tidal constituents could also play a role. It can be noticed that such oscillation can happen due to tide surge interaction phenomenon which tends to modify the tidal phase (Flather 2001). Such phenomenon happens in the presence of significant tide range (and induced tidal currents) and surge. Saint-Malo is characterized by a strong tidal range (about 10.7 m for a mean spring tide) and annual surges of 0.8 m. This is thus a site which can be subject to this type of phenomena (Idier et al. 2012a). This can also be due to harmonic analysis made on old water level measurements (before numeric tide gauges). There are sometimes small phase shift errors on measurements due to mechanic system, generating phase shift in harmonic components. The consequence can be this kind of oscillations showing a phase shift between measured and predicted water levels. For these reasons, in our study, we have preferred to use only data from numeric gauges. Another explanation could be a small harmonic constituent which has not been taken into account into the harmonic analysis. This is quite unlikely, because the harmonic constituent list is quite complete (143 constituents). All the same, such oscillations are not systematic, and the same kind of graph at Dunkerque does not show them (not presented here).

Results
Three different meteorological forcing have been tested (Fig. 20): & ARPEGE High Resolution merged with AROME, which corresponds to spatial interpolation between ARPEGE HR and AROME over the domain & ARPEGE High Resolution & ARPEGE The merge consists of a spatial interpolation of zonal and meridional wind components from these datasets on a common grid, whose spatial resolution is the smallest of the two model resolutions (0.025°for AROME). The spatial interpolation is bicubic. A relaxation mechanism is applied, in order to allow a progressive transition between the two datasets. This transition is linear and is applied on a 2°zone around the high spatial resolution model (AROME).
Modeled surges are compared with measurements and also with surges from current operational model, labeled as "old operational model" on figures, which runs in PREVIMER since 2006. This operational model is similar to the reference configuration (with a constant drag coefficient equal to 0.0016, and ARPEGE meteorological forcing), but its spatial resolution is of 5.6 km instead of 2 km.
The influence of meteorological forcing on surges in December 2011 and January 2012 is presented on Figs. 21 and 22. The highest storm surge over this period reaches almost 2 m at Dunkerque on the 5th of January 2012. Surges from the 4th to 7th of January 2012 are presented on Fig. 23. At Saint-Malo, measured surges present typical oscillations, which are discussed in Sect. 4.2 (Fig. 19). For the different meteorological forcing, statistics have been computed for January 2012: -Root mean square error (RMSE) values for surges are presented in Table 5. -Maximum surges are presented in Table 6.  -Error indicators, mentioned in Sect. 3, have been computed for January 2012 and for Andrea storm from the 4th to the 7th of January. Indicators are very similar for these two periods; therefore, only indicators computed from the 4th to the 7th of January are presented (Fig. 23).
First of all, these results show that there is no model improvement in using ARPEGE HR merged with AROME rather than ARPEGE HR alone. The results are very similar with for example at Dunkerque RMSE equals to 12 cm for ARPEGE HR merged with AROME and 11 cm for ARPEGE HR. Moreover, the RMSE is generally smaller for ARPEGE HR alone, but AROME improves slightly the peak error. AROME has the same time resolution as ARPEGE HR but a better space resolution; sensitivity tests carried out on wind forcing in the Gulf of Lions (Mediterranean Sea) suggested the prevailing of space resolution on time resolution (Schaeffer et al. 2011).
Furthermore, RMSEs for ARPEGE HR and ARPEGE are also very similar (respectively 11 and 10 cm at Dunkerque), Fig. 18 a Atmospheric pressure, b wind direction, and c wind speed measured (data) and modeled (AROME, ARPEGE HR, and ARPEGE meteorological models) from the 4th to 8th of January 2012 at Dunkerque Comparing with the old operational model, even with the same meteorological forcing (ARPEGE), the maximum surge reaches 1.45 m instead of 1.27 m. This is not only due to a better spatial resolution (2 km instead of 5.6 km) but also due to a better sea surface drag parametrization taking into account wind and wave with a Charnock parameter derived from WAVEWATCH III® model. However, even improving the meteorological forcing and the sea surface drag parameterization, the model for this event (5th of January, at Dunkerque) still tends to underestimate maximum surge (1.62 instead of 1.93 m) (see Fig. 23). This could be explained by wave setup at tide gauge, which is not modeled here. Thanks to modeling, Bertin et al. (2012b) estimate a wave setup of 40 cm during Xynthia storm where the coastline was directly exposed to the storm. There is also a phase shift between model and measurements, which could be explained with input data (bathymetry), parameterization (bottom friction is here constant all over the domain), and forcing (tidal model at boundary and meteorological forcing). Moreover, model resolution cannot reflect the details of the bathymetry of the studied area (for instance, in Dunkerque, there are sand banks which evolve with time). The wave setup can also impact the phase shift since its maximum is not necessarily reached at the same time as the atmospheric surge.

Conclusion
This paper aimed at better understanding and reproducing the storm surge dynamics along the French and English Channel coasts, focusing on (1) sea surface drag parameterization and (2) meteorological forcing data quality. Within this framework, two sensitivity studies have been carried out. The first one showed some limits of only wind-dependent sea surface drag parameterizations for storm surge modeling. Indeed, four wind-dependent parameterizations (Moon and Makin formulations, Charnock formulation with a constant Charnock parameter (α c =0.014 and α c =0.022), and one wind-and wavedependent parameterization (Charnock formulation with a variable α c )) have been tested. Storm surges issued from modeling based on these five formulations were compared to those computed with a reference configuration that uses a constant drag coefficient (0.0016) and also to observations. The comparative study has revealed the added value of using the Charnock formulation for storm surge modeling. The Charnock formulation taking into account wave conditions (variable α c ) enabled strongly reducing errors, by reproducing higher surge peaks compared to the other formulations, with increase from several centimeters to dozen of centimeters for each studied locations (Dunkerque, Saint-Malo, Le Conquet, and La Rochelle). It should be noticed that the use of a high constant α c gave similar results but with possible side effects. However, the limited impact of a variable α c compared to a constant one set to 0.022, in our case, is due to the particular parameterization of the wave model. As a result, the use of a   variable α c rather than a constant one depends on its link with wave age and therefore on the parameterization used in the wave model to compute it. A variable α c that translates the presence of steep and young wind-waves with more realism should induce amplified storm surges. The second sensitivity analysis puts the stress on the importance of the meteorological forcing quality used for storm surge modeling. Meteorological forcings with different time and space resolution have therefore been tested. Results showed that ARPEGE High Resolution improved clearly maximum surges, up to 25 cm, compared with ARPEGE forcing. Thus, improvements of meteorological forcing time and space resolution, parameterization of drag coefficient taking into account wind and wave, and also the space resolution of the ocean model itself allowed to divide by two water levels RMSE (40 cm to 22 cm) and to improve maximum surges estimation up to 35 cm (127 to 162 cm at Dunkerque during Andrea storm). Nevertheless, the model still underestimates maximum surges; for example, measured surge reaches 193 cm at Dunkerque for Andrea storm. At this stage, it is important to underline that model outputs and measures are not rigorously comparable because of tidal gauge locations (in harbors) and not modeled physical phenomena such as wave setup and seiches. However, seiches can be filtered out in the tidal gauge measurements (this was the case in the sensitivity analysis to the drag formulation). For the wave setup, at harbor locations, its magnitude is much smaller than on surrounding coasts where it can reach several tenths of centimeters on open beaches (Idier et al. 2010). For instance, at La Rochelle and for the Xynthia event, Bertin et al. (2012b) provide an estimation of 5 cm. Such value falls in the error range obtained between the modeled and observed storm surge. Thus, within the present paper, we reduced the storm surge underestimation of observation by the model and obtained a so-called best model providing an underestimation of the storm surge consistent with the wave setup magnitude order. For further refinement of the storm surge model, it would be worthwhile to properly estimate, probably by modeling, the wave setup at the harbor scale. Moreover, this underestimation could be partially explained by the quality and resolution of bathymetry (particularly in areas where bedforms-dunes and sandbanks-are present), the bottom drag friction parameterization which does not take into account wave action and wave setup. In addition, the 3D dynamics which is not taken into account could also contribute to explain this underestimation. Indeed, a basic difference between 3D and 2D models for storm surge modeling is the bottom friction parameterization (for more details, see Weisberg and Zheng 2008). As shown in many applications (e.g., Davies 1988;Weisberg and Zheng 2008), the bottom stress is generally overestimated in the 2D case, such that the pressure gradient force tends to be underestimated. Hence, the storm surge heights are also underestimated by the 2D model relative to the 3D model. Moreover, the quality of the Charnock parameter used in the sea surface drag formulation also contributes to the modeled storm surge quality. Indeed, if the wave model underestimates the sea surface roughness, then the Charnock parameter is also underestimated and the modeled surge is weaker (Bertin et al. 2012a).