1 Introduction

Storm surges can have enormous consequences, like flooding of areas adjacent to seas and oceans. Buildings and nature are at risk of heavy damage, and societies can be seriously disrupted. At mid-latitudes, for example around the North Sea, coasts have always been subject to the risk of flooding due to surges generated by heavy storms. The Netherlands, Great Britain and Germany have experienced several of these severe floods throughout history. At times, it caused thousands of casualties (Brooks and Glasspoole 1928). The most recent big flooding in the south-western part of the Netherlands occurred on 1 February 1953.

Tropical cyclones can also cause devastating storm surges, a well known example being Hurricane Katrina in 2005 where, besides the damage that was caused by strong winds, about 80% of the city of New Orleans was flooded due to the accompanying surge. Other examples of regions that are vulnerable for storm surges induced by tropical cyclones are the coasts on the Bay of Bengal and the Philippines.

In the second half of the twentieth century, the development of numerical models facilitated the prediction and investigation of storm surges. Most of such models use non-linear depth-averaged shallow water equations (cf. WMO 2011, and references therein). Due to the introduction of these models, knowledge about storm surges and of the mechanisms behind them has increased enormously. This, together with the development of atmospheric circulation models, resulted in relatively accurate real-time forecasts of storm surges, which help to prepare for such events.

A major error source in the prediction of storm surges is the uncertainty in the meteorological input. Storm surge models are forced by wind and pressure fields, produced by atmospheric models. Errors in these input fields are transferred into errors in surge height along the coast. An aspect that has not been studied much in this context is the influence of the short-term variability of the wind. Zeng et al. (2002) derive a parametrisation for wind gustiness for the computation of ocean surface fluxes, especially for the tropics, and conclude that it is important mainly under weak wind conditions. Koopmann (1962), however, suggests that the surge level is enhanced in cases with very high gustiness. Cavaleri and Burgers (1992) investigated the consequences of wind gustiness on the growth of wind waves. They found that growth rate and energy of the waves are affected by short-term wind variability. Wave energy was increased for high values of gustiness. Analogous to its effect on waves, we expect an effect of short-term wind fluctuations on storm surges, because τ ∼ u 2, where τ is the instantaneous wind stress and u the wind velocity. Thus, fluctuations in the wind speed contribute to the mean stress.

In this paper, we investigate whether part of the error in storm surge model forecasts for the North Sea can be explained by the effect of short-term wind fluctuations and assess the importance of this contribution compared to other error sources, for example in atmospheric forcing or bathymetry.

The storm surge model that has been used for the experiments is the Dutch Continental Shelf Model (WAQUA/DCSM). It solves the depth-averaged non-linear shallow water equations on the NW European Continental Shelf, and it is based on the WAQUA-in-SIMONA software package, supplied by Rijkswaterstaat (Gerritsen et al. 1995; Verlaan et al. 2005). The model is used for day-to-day forecasts of water levels along the coast of The Netherlands with input from KNMI’s Hirlam model (HIRLAM Consortium 2009) and the ECMWF-IFS model (ECMWF 2004).

As a measure for the amount of short-term wind fluctuations, we will use gustiness, which we introduce in Sect. 2 and define as the standard deviation of the (longitudinal) short-term fluctuations. This term should not be confused with the related wind gust, which is the maximum of the wind in a prescribed time interval. This is an output parameter of atmospheric models like ECMWF-IFS or Hirlam. In this study, we will deduce the gustiness from the wind gust.

To investigate the exact effect of varying gustiness on the storm surge, we designed experiments with simple idealised wind fields. In Sect. 3.1, random terms are added to slowly varying uniform winds to model short-term wind fluctuations, and sensitivity experiments are performed. In this way, the effect of gustiness on storm surge characteristics can be studied precisely in controlled situations.

To incorporate the effect of short-term wind variations into a storm surge model, one has to transform the stochastic gust output as produced by an atmospheric model into a deterministic modification of the wind stress. From the results of the idealised sensitivity experiments, we extract a parametrisation for varying gustiness, which is presented in Sect. 3.2. With this parametrisation, we performed some case studies of severe storms in the North Sea in 2007 and 2008, using wind, pressure and wind gust from the ECMWF-IFS model as input.

Finally, in Sect. 4, we present and analyse a long run for the entire storm season September 2007 to April 2008. For four stations along the Dutch coast: Vlissingen, Hoek van Holland, Den Helder, and Harlingen, we investigate the effect of the inclusion of varying gustiness on the errors in the forecasts of high and low tides.

One has to keep in mind that the mean effect of the short-term variations is already taken into account in storm surge models. This is due to the tuning of the drag relation, when the model is calibrated for a certain domain to optimise the correspondence of the observations and model outcomes. In this way, a mean contribution due to short-term variations is already in the model. Nevertheless, from situation to situation, the amount of the short-term fluctuations, and hence their contribution to the stress on the sea surface, may vary. This means that adding gustiness explicitly to the model forcing will overestimate its net effect on the surge levels and would initially deteriorate the model results. Therefore, the drag relation would need to be retuned to remove the effect of the average gustiness.

2 Short-term fluctuations of the wind

Traditionally, operational storm surge models are forced with, for example, hourly wind fields created by atmospheric models. These wind fields are interpolated in time, as depicted by the red line in Fig. 1. However, the figure also shows a time series of observed wind speed, which varies on much smaller time scales, and those short-term fluctuations are not accounted for in the model. Short-term motions in the lowest part of the atmosphere occur in complicated irregular patterns, both in space and time, but can be described in statistical terms, as is explained in, for example, Tennekes and Lumley (1972) and Garratt (1992).

Fig. 1
figure 1

Time series of observed wind speed (green, fluctuating) and the corresponding hourly input for a storm surge model (red)

Full size image

The stress that is exerted on the sea surface by the wind can be expressed as a quadratic function of the horizontal wind, the drag relation, and also points in the direction of the wind:

$$ {\boldsymbol{\tau}} = \rho_{\rm a}\;C_{\rm d} \; \|{\bf u}\| {\bf u}. $$
(1)

Here, \({\boldsymbol{\tau}}\) is the stress vector, \(\rho_{\rm a}\) is the air density, and u is the wind vector at 10 m height. The drag coefficient C d is a bulk parameter, which depends on atmospheric conditions and the wave state. For neutral conditions, it is described by the Charnock relation (Charnock 1955). As the wind speed enters quadratically into the stress, and even also into the drag coefficient, short-term fluctuations around the mean value do not average out. Thus, the average stress is altered by variations in the wind speed on small time scales. In this study, we differentiate between slowly varying wind \(\overline{{\bf u}},\) with time scales of one or more hours and longer, and faster variations u′. The total wind vector then is \({\bf u} = \overline{{\bf u}} + {{\bf u}}'.\)

The slowly varying wind corresponds to the synoptic scale wind patterns that are provided by atmospheric models and serve as traditional input for the storm surge model. The time scale depends on the availability of model output, which is, for example, 3 h for the ECMWF-IFS model and 1 h for Hirlam. The mean effect on the sea level of variations with time scales shorter than the input time step is corrected for via tuning of the model. But deviation in space and time from the mean is not taken into account and can be considered unresolved motion.

To assess the contribution of variations on these smaller time scales, we express the vectors as a longitudinal component in the direction of \(\overline{{\bf u}}, \) and a lateral component:

$$ \overline{{\bf u}} = {\underline{\underline{{\mathbf R}_\phi}}} \; \left( \begin{array}{l} U \\ 0 \end{array} \right), \quad \quad {\bf u}' = {\underline{\underline{{\mathbf R}_\phi}}} \; \left( \begin{array}{l} u' \\ v' \end{array} \right). $$
(2)

Here, u′ and v′ are longitudinal and lateral components of short-term wind fluctuations, U is the mean wind speed, and

$$ {\underline{\underline{{\mathbf R}_\phi}}} = \left( \begin{array}{ll} -\sin\phi& \cos\phi \\ -\cos\phi & -\sin\phi \end{array} \right) $$

is the rotation tensor for the mean wind direction ϕ, defined in Fig. 2.

Fig. 2
figure 2

Coordinate system with total wind u (black) split up in mean wind \(\overline{{\bf u}}\) (blue) and short-term deviations u′ (red). The wind direction ϕ is defined as ‘coming from’ with respect to North

Full size image

The time averages over the short-term fluctuations are

$$ \overline{u'} = 0, \quad {\overline{u^{'2}}}= \sigma_{u'}^{2}, \quad \overline{v'} = 0, \quad {\overline{v^{'2}}} = \sigma_{v'}^{2}, \quad \overline{u'v'} = {\mathop{\rm{cov}}}(u',v'), $$
(3)

where σ u′ and σ v′ are standard deviations of short-term longitudinal and lateral fluctuations and \({\mathop{\rm{cov}}}(u',v')\) is the covariance between u′ and v′.

To calculate the mean stress that is exerted on the sea surface, Eq. 2 is substituted into Eq. 1 and averaged in time. With a second-order Taylor expansion around u′ = 0 and v′ = 0, and the assumption that O(u′) ≈ O(v′) ≪ O(U) and that O(u′)3 and O(v′)3 are small (Smith and Chandler 1987), the resulting mean stress vector reads

$$ \begin{aligned} \overline{\boldsymbol{\tau}} &= \overline{\rho_{\rm a}\;C_{\rm d}\; \sqrt{\left (U + u'\right )^2 + {v'}^2} \;\; \underline{\underline{{\mathbf R}_\phi}} \; \left(\begin{array}{l} U + u' \\ v' \end{array}\right)} \\ & \simeq \rho_{\rm a}\;C_{\rm d} \; \underline{\underline{{\mathbf R}_\phi}} \; \left(\begin{array}{l} \overline{U^{2}} + \overline{u^{'2}} + \overline{2Uu'} + \overline{{1 \over 2} v^{'2}} \\ \overline{Uv'} + \overline{u'v'} \end{array}\right) \\ & = \rho_{\rm a}\;C_{\rm d}\; \underline{\underline{{\mathbf R}_\phi}} \; \left(\begin{array}{l} U^2 + \sigma_{u'}^2 + \frac{1}{2}\sigma_{v'}^2 \\ \mathop{\rm{cov}}(u',v') \end{array}\right). \end{aligned} $$
(4)

Hence, the components of the stress are enhanced with terms that contain the standard deviations and covariance of the components of the short-term wind fluctuations. Note that also the drag coefficient depends on the wind speed. It increases for higher winds and thus increases the stress even more. Additional terms that result from this dependence are not accounted for in Eq. 4.

In order to predict the value of the increased wind stress, we need to derive or specify expressions for the unknowns in Eq. 4: the second-order moments σ u′,  σ v′, and \(\mathop{\rm{cov}}(u',v'). \) We assume that the correlation between the short-term wind components is small, \(\mathop{\rm{cov}}(u',v') \approx 0\). As a result, the stress remains pointing in the original wind direction. Furthermore, the standard deviations of the short-term fluctuations can be linked to each other:

$$ \sigma_{v'} = r \; \sigma_{u'}. $$
(5)

The ratio r is a constant, as follows from considering density spectra of the variation velocity components. Kaimal et al. (1972) deduce from the 1968 AFCRL Kansas experiments that under neutral conditions

$$ \frac{\sigma_{v'}}{\sigma_{u'}} \approx 0.75. $$
(6)

Now, only one unknown is left to predict the extra stress due to short-term deviations. This unknown σ u′ can be extracted from atmospheric models, such as ECMWF’s IFS (see Sect. 4). In this study, we will thus take the standard deviation of the longitudinal wind fluctuations as a measure for short-term variations in the atmosphere, and call it gustiness.

As already mentioned in the Introduction, due to the tuning of the drag relation, the average effect of short-term motion is already accounted for. What is not properly accounted for, however, is variation in the short-term motion, that is, higher or lower σ u′ in space and time.

3 Experiments

The impact of the enhancement of the wind stress by gustiness has been studied with the WAQUA/DCSM model in a series of Monte Carlo experiments to determine the sensitivity of the surge level to different values of σ u′.

For this purpose, we take a mean wind forcing that resembles the forcing that is used by Gill (1982, pp. 394–403), where he gives an analytical solution of the February 1953 storm in the North Sea. A uniform north–westerly wind field is created, parallel to the orientation of the North Sea, that is, a wind direction of ϕ = 315°. This is the most effective direction for high storm surges in the Netherlands, because of the long wind fetch. The mean wind speed U of this wind field varies in time, according to a harmonic sine function with a period of 6 days. After 3 days, the wind speed is set to zero, thus

$$ U({ \hbox {\bf x}}, t) = \left\{ \begin{array}{lll} M \sin(\frac{t \pi}{t_s})& \hbox{when} &\; \; 0 \leq t \leq t_{\rm s} \;, \\ 0 & \hbox{when} &\; \; t > t_{\rm s} \;, \end{array} \right. $$
(7)

with t time, M maximum mean wind speed, here M = 25 m/s and t s is duration of the storm: t s = 3 days. The initial state of the sea level is chosen to be the astronomical tide on 1 January 2009 0 UTC. This state is achieved by running the model over the period of a month without meteorological forcing.

3.1 Monte Carlo experiments

Gustiness is added to the wind forcing in a stochastic way. We assume that the longitudinal and lateral velocity components of short-term fluctuations are Gaussian distributed (as in a neutral boundary layer, Chu et al. 1996). This corresponds to the assumption that higher than second-order moments are zero in Eq. 4. We define \(\vec{u}_s\)(x, t) on grid points \({\bf x} = (\lambda_i, \varphi_j), \) as an N-dimensional state vector \(\vec{u}_s\)((j + Bi), t). Here, N is the total number of grid points, λ i and \(\varphi_j\) are longitude and latitude of the grid point, and i and j are column and row of the grid point, counted from the south--west corner. The columns of \(\vec{u}_s\) are placed below each other in the state vector, with B the number of grid points along the meridional direction. Now, we model

$$ \vec{u}_s = \underline{\underline{{\mathbf D}_u}} \cdot {\vec{\eta}_u}, $$
(8)

where \(\underline{\underline{{\mathbf D}_u}}\) is an N × N tensor and \({{\vec{\eta}}_{u}}\) is an N-dimensional vector. The components of \({{\vec{\eta}}_{ u}}\) are drawn from a random number generator and have a Gaussian distribution with mean 0 and standard deviation 1. With \(\underline{\underline{{\mathbf D}_u}}, \) the standard deviation of the gustiness can be manipulated and spatial correlation can be introduced. The same procedure is used to model gustiness in the lateral direction \(\vec{v}_s\).

With U(x, t) from Eq. 7, the total wind forcing in grid point x and at time t, is then

$$ {\bf u}\,({\bf x}, t) = \underline{\underline{{\mathbf R}_\phi}} \; \left( \begin{array}{l} U({\bf x}, t) + \vec{u}_s({\bf x}, t)\\ \vec{v}_s({\bf x}, t) \end{array} \right). $$
(9)

Figure 3 shows a typical wind speed at an arbitrary grid point in a uniform wind field, created with Eq. 9.

Fig. 3
figure 3

Example of wind speed with no added gustiness (red, smooth) and random gustiness with σ u  = 0.2U (blue)

Full size image

3.1.1 Gustiness as random noise

From Eq. 8, we generate uncorrelated short-term wind fluctuations by the addition of random noise. To obtain noise with standard deviations σ u′ and σ v′, we set

$$ \begin{aligned} \underline{\underline{{\mathbf D}_u}} &= \sigma_{u'} \underline{\underline{\mathbf I}} \quad \hbox{and}\\ \underline{\underline{{\mathbf D}_v}} & = \sigma_{v'} \underline{\underline{\mathbf I}}, \\ \end{aligned} $$
(10)

where \(\underline{\underline{\mathbf I}}\) is the identity tensor. We take gustiness proportional to the mean wind speed, that is, σ u′ = α u U(t) and σ v′ = α v U(t). According to Eq. 6, we take the ratio of the standard deviation in the lateral and longitudinal direction 0.75 and choose α u  = 0.4 and α v  = 0.3. These are relatively high values for the standard deviation of short-term fluctuations.

With these settings, we create an ensemble of 50 realisations to examine the sensitivity of the sea level to gustiness. For every model time step and for each of the realisations, a new set of random numbers is drawn for \({{\vec{\eta}}_{u}}\) and \({{\vec{\eta}}_{v}},\) and the storm surge model is run with Eq. 9 as input. Figure 4 gives an example of a time series from the ensemble for Hoek van Holland, one of the important locations along the coast of The Netherlands. The red line gives the surge without gustiness. The blue and orange lines are the 50 members of the ensemble. Clearly, the surge is amplified by the addition of gustiness.

Fig. 4
figure 4

Surge in Hoek van Holland of 50 ensemble runs with random gustiness in longitudinal and lateral direction: α u  = 0.4 and α v  = 0.3 (orange lines and blue line). Red line, separated from the others, is the control run without gustiness

Full size image

3.1.2 Gustiness as correlated noise

Although the experiment from Sect. 3.1.1 gives an qualitative insight in the role of gustiness in storm surges, it is far from realistic. We further investigated the influence of gustiness by adding correlation in both space and time. Correlation in space is created by introducing non-zero off-diagonal elements in the matrix \(\underline{\underline{{\mathbf D}_u}}\). We use a weighted average over a zone of 80 km around each grid point with weights decreasing linearly with distance and renormalised to maintain the same variance in the gustiness field as before.

Correlation in time is created by shifting the gustiness pattern and adding it to the next time step:

$$ \vec{u}_s(\lambda_i, \varphi_j, t_k) = \varepsilon \;\; \vec{u}_s(\lambda_{i'}, \varphi_{j'}, t_{k-1}) + (1 - \varepsilon) \; \underline{\underline{{\mathbf D}_u}} \cdot {{\vec{\eta}}_{u}}(t_k), $$
(11)

where the translation \((\lambda_{i'}, \varphi_{j'}, t_{k-1}) \rightarrow (\lambda_i, \varphi_j, t_k)\) corresponds to a large scale movement of small-scale weather patterns, like, for example, showers. We used \(\varepsilon = 0.75\) and a translation of 30 m/s parallel to the average wind direction to mimic a typical system that moves with 100 km/h.

The overall conclusion from these experiments is that the average of the storm surges from the ensembles is not sensitive to adding correlation in the gustiness. In contrast, the width of the ensemble band increases, which seems admissible as the effect of correlated gustiness in neighbouring grid points and time steps would accumulate whereas uncorrelated gustiness would partly neutralise.

3.1.3 Short-term wind fluctuations with higher cut-off frequencies

We also experimented with the time scales of the added wind variations by reducing the model time step. If gustiness is fed into the model with smaller time intervals, the cut-off frequency of the short-term fluctuations is increased.

We compared ensembles created with random gustiness with time steps of 10 min (the standard in WAQUA/DCSM) and 1 min, but identical σ u′ and σ v′. Both ensembles did not differ significantly, and we concluded that the higher frequencies in the gustiness do not enhance the surge in WAQUA/DCSM.

3.2 Parametrisation

Up to now, we have used a slowly varying wind field U that is uniform in space, as is the standard deviation of the short-term fluctuations. In realistic situations, this is not the case, as wind and gustiness will vary in space. Also, we want to eliminate the ensembles for practical applications. This prompts for a parametrisation of the gustiness effect on storm surges.

As suggested by Eq. 4, we propose to replace the wind forcing from Eq. 9 by

$$ {\bf u}({\bf x},t) = \underline{\underline{{\mathbf R}_\phi}} \; \left( \begin{array}{c} \sqrt{U({\bf x},t)^2 + \beta \sigma_{u'}({\bf x},t)^2} \\ 0 \end{array} \right). $$
(12)

The longitudinal wind U(x, t) + u′(x, t) is replaced by the quadratic addition of the mean wind speed and a factor β times the standard deviation of the longitudinal short-term fluctuations. Corresponding to the assumption that \(\mathop{\rm{cov}}(u',v') \approx 0\) in Eq. 4, the lateral component of the wind in this parametrisation is 0. This representation is convenient for practical applications later on, because both U and σ u′ are output parameters of atmospheric models that vary in space and time. Furthermore, β is a constant. We use it to tune the results with the new wind forcing to the results with the ensembles.

The longitudinal component in Eq. 12 can be rewritten as

$$ \sqrt{U({\bf x},t)^2 + \beta \sigma_{u'}({\bf x},t)^2} = U({\bf x},t) \sqrt{1 + \beta \alpha_u^2}, $$
(13)

where α u  = σ u′/U is the ratio of the short-term variations and the slowly varying wind speed, as in Sect. 3.1.1. Note that in this parametrisation only the longitudinal standard deviation is used. As already discussed in Sect. 2, the lateral standard deviation is assumed to be a factor r times the longitudinal standard deviation, that is, α v  = rα u , with r = 0.75 (Kaimal et al. 1972). We will compare the outcomes of the experiment in Sect. 3.1.1 (α u  = 0.4, α v  = 0.3) with the results of the run of the model with the new parametrisation for the wind forcing. We thus use α u  = 0.4. We do not include correlation in the comparison, because the average of the different realisations does not deviate substantially from the runs with random noise.

Figure 5 summarises the validation of the parametrisation for this configuration. The value of β was tuned to 2.2 to obtain the same ensemble-mean surge height with Eq. 12 as that obtained with added random gustiness (see Fig. 4).

Fig. 5
figure 5

Surge in Hoek van Holland of 50 ensemble runs (orange) and the parametrisation (black), with β = 2.2. Red line, separated from the others, is without gustiness

Full size image

To get more confidence that the parametrisation properly describes the effect of gustiness on a storm surge, we also conducted a run with a more complicated wind field. In this experiment, we divided the North Sea in 4 sectors, each with their own values for M, ϕ (Eq. 7) and α u . In the NW sector, the wind comes from the Northwest with M = 30.0 m/s and α u  = 0.32, in the NE sector from the Northeast with M = 37.5 m/s and α u  = 0.16, in the SW sector the wind is South with M = 12.5 m/s and α u  = 0.40, and in the SE sector the wind comes from the Southwest with M = 10.0 m/s and α u  = 0.60. Also in this experiment, the surge with the parametrised gust with β = 2.2 lies within the ensemble band, although not exactly in the centre anymore.

4 Application

For the application of the parametrised gustiness in real situations, we use meteorological input from the ECMWF model archive. One of the output parameters of this model is the wind gust, the maximum 3-s wind speed to occur within a 3-h interval at a height of 10 m (see Beljaars et al. 2001).

To translate the wind gust into the gustiness as used in the previous sections, we use the probability P that the measured wind gust in a time period T exceeds a value G. This is a function of the distribution of the wind, characterised by the mean wind speed U and longitudinal standard deviation σ u′ (gustiness), and the time scale of short-term deviations τ. Characteristic time scales for short-term deviations can be found again from the spectra of Kaimal et al. (1972). When imposing a value for P,  σ u′ can be calculated from G = U + gσ u′ (Beljaars 1987; Wichers Schreur and Geertsema 2006), with

$$ g = \left(2\ln\left(\frac{T}{\tau}\frac{1}{\sqrt{2 \pi}\ln(\frac{1}{1-P})}\right)\right)^\frac{1}{2}. $$
(14)

The resulting value G can be interpreted as a prediction for the wind gust, with an exceedance probability P. In the ECMWF-IFS model, an exceedance probability of 25% and a time period T of 3 h is used, which lead to g = 2.93.

Thus, in our runs, of which the results are presented below, we use G from the ECMWF-IFS model to compute σ u′.

4.1 Individual cases

The influence of the inclusion of parametrised gustiness on the calculated storm surge was investigated for a number of individual storm cases. We selected a large surge on 9 November 2007, a surge on 7 December 2007, and a negative surge in January 2008. Figure 6 gives the time series of the total sea level during these events in Hoek van Holland.

Fig. 6
figure 6

Time series of the sea level in Hoek van Holland for selected cases: Astronomical tide (dotted blue), observed (green dots), calculated without (solid red) and with (dashed black) explicit gustiness

Full size image

Generally, gustiness increases the surges, but the extent varies. In Fig. 6, on the first peak on 9 November, gustiness improves the model result, but otherwise, where the modelled sea level deviates substantially from the observed sea level, the inclusion of gustiness does not close this gap.

4.2 Analysis of a winter season

For a statistical analysis of the impact of gustiness on calculated storm surges, the storm surge season 1 September 2007–1 April 2008 was selected. During this winter, several storm surges occurred, and also a period of moderate negative surges.

As the average effect of gustiness is included in the tuning of the air–sea interaction, no overall direct improvement of the forecasts should be expected from adding it explicitly. Instead, adding gustiness will always increase the calculated surge, less so in cases where gustiness is low and more when it is high, but the model was originally tuned to be unbiased on average.

We expect that the original model underestimates the surge when gustiness is high and overestimates with low gustiness. Hence, we will be looking for a relation between the forecast error \(\Updelta H,\) the difference between the calculated surge H mod and the observed surge H obs, and the added gustiness. A correlation between gustiness and error would imply that gustiness explains part of the error in the storm surge forecasts. If we define \(\Updelta H = H_{\rm mod} - H_{\rm obs},\) this correlation should then be negative.

We investigate the possible presence of a relationship between \(\Updelta H\) and gustiness by examining high and low tide skew surges (see Gerritsen et al. 1995). The gustiness effect is quantified by the difference between the calculations with and without added gustiness, and the error \(\Updelta H\) is taken from the calculation without gustiness. As an example, Fig. 7 shows the results for four gauges along the coast of The Netherlands: Vlissingen, Hoek van Holland, Den Helder, and Harlingen.

Fig. 7
figure 7

Forecast error \(\Updelta H\) versus gustiness effect for high and low tide skew surges and positive (red pluses) and negative (blue squares) surges

Full size image

The figures do not show statistically significant correlations. The weak signals for, for example, high tides in Vlissingen or low tides in Den Helder, rather mirror a similar correlation between the error and the skew surge itself (not shown), which in turn is correlated to the gustiness effect. For other locations where the astronomical tide is less accurately known, this can be even bigger and would completely overwhelm any effect of the gustiness itself.

Another observation from Fig. 7 is that the error in the forecasts is large compared to the error that is solely caused by the effect of gustiness. The range of the forecast error is in the order of 50 cm, and the maximum addition due to gustiness is only 10 cm. This means that before the role of gustiness could become important in storm surge forecasting, errors due to other effects, for example in the meteorological forcing or bathymetry that is used in the model, have to be decreased first.

5 Discussion and conclusions

We investigated the influence of gustiness on storm surge calculations in the North Sea. With Monte Carlo experiments, we examined the sensitivity of the sea level to gustiness and used this to derive a parametrisation for calculations with real-world input.

In hindcasts of a few individual cases and a complete winter season, we did not find a clear sign that the explicit inclusion of gustiness helps to decrease the differences with observations. Moreover, the error that is solely caused by the effect of gustiness is small compared with the overall forecast error.

To investigate the sensitivity of the results to the choices, we made in the parametrisation, we repeated the hindcast of the winter season with strongly increased gustiness, by setting β = 5 instead of β = 2.2, as follows from Sect. 3.2. The resulting plots of \(\Updelta H\) versus gustiness effect are almost identical to those from the previous experiment, with a scaled gustiness axis. This means that the exact choice of the parameter β does not affect the conclusion that gustiness does not explain the error in storm surge calculations. Gustiness could be important in other domains, but for the North Sea the effect on surges is far too diffuse to find an improvement.

Hence, the main conclusion of the study is that, at the moment for the North Sea, it is not crucial to incorporate the effect of gustiness explicitly in storm surge models. However, as the error in the calculations will be reduced in the future, for example by improvement of the meteorological forcing, it can be useful at some moment to investigate the effect of the addition of explicit gustiness on surges again.