Storm surge modelling

Storm surge models are generally based on the two-dimensional shallow water equation to compute the water level and velocity. According to different modelling purposes, the storm surge model can be divided into: i) storm surge model without tide; ii) storm surge model including astronomical tide; and iii) storm surge flooding model considering inundation.

In a storm surge model without tide, only the effect of the meteorological forcing needs to be considered. Obtained result is only the water level rising and falling under the effect of atmospheric pressure and wind. The role of this model is generally to provide forecasters with a reference for the magnitude of storm surge when making forecasts, and it cannot truly reflect the fluctuations of water level.

Based on the former model, the tidal boundary conditions and tidal potential are introduced in the total water level storm surge model, in which the nonlinear interaction between storm surge and astronomical tide can be fully considered in the model. This model can be used for releasing storm surge numerical forecast products and providing reference for coastal response to storm surge disasters. It can also provide the water level changes in the target area under the influence of extreme weather conditions for coastal engineering (such as harbours, wharfs, seawalls, offshore wind farms, etc.).

The storm surge flooding model considering inundation is more complicated than the previous two models. The interaction of storm surge and astronomical tide, and the interaction of storm surge and wave are also considered in the model. The inundation range and depth of the coastal area can be obtained by a storm surge simulation. Figure 7.9 shows the main storm surge modelling components used by a forecasting system and that will be detailed in the next sections.

Reliable and accurate shoreline and bathymetry data are the basis for modelling of storm surge, tidal, and storm surge flooding models. In the process of model setup, the computational grid should be determined according to the demand, and then the data of shoreline and bathymetry should be collected according to the location and scope. 

Sources of bathymetric data can be found on Section 4.2.4. These data need to be used with caution when establishing high-resolution storm surge models. As the bathymetry of coastal or estuary changes rapidly over time, these upstream data may not be able to be updated in time. Therefore, the correctness of the data needs to be verified before it is used.

When the published upstream data do not meet the requirements, it is preferable to use the latest sea chart or high-resolution DEM data. In addition, the datum of the data needs to be confirmed as, due to different sources of data, the datum could be different. In order to truly reflect the effects of bathymetry on storm surges, data from different sources need to be unified on the same datum.

In the storm surge and storm surge flooding modelling, tidal waves are generally used as open boundary conditions. Open boundary conditions can be velocity, water level, or harmonic constant. 

Tidal height at any location and time can be written as a function of harmonic constituents according to the following general relationship:

Where: 

• H(t) = height of the tide at any time t
• H_o = mean water level above some defined datum such as mean sea level
• H_n = mean amplitude of tidal constituent n
• f_n = factor for adjusting mean amplitude H_n values for specific times
• a_n = speed of constituent n in degrees/unit time
• t = time measured from some initial epoch or time, i.e., t=0 at t_o
• (V_o+ u)_n = value of the equilibrium argument for constituent n at Greenwich and when t=0
• g_n= epoch of constituent n, i.e., phase shift from tide-producing force to high tide from t_o

In Eq. 7.1, for a tidal component, the f_n and the (V_o+ u)_n will change with time, and the H_n and the g_n will change with the geographical location. Therefore, according to the start time of simulation, the fn and (V0+ u)n of the proposed tidal component should be set in the model, and the Hn and the gn boundary conditions should be given according to the position of the boundary grid node.

The location of the tidal open boundary is critical. First of all, it is necessary to ensure that the grid can completely cover the definition area. Secondly, it is better to set the open boundary near the tidal station, because the tidal constituent data near the tidal station is more accurate. Finally, it is preferable not to set the open boundary at the amphidromic point nearby, because too small amplitude will bring uncertainty to the simulation.

The data for tidal boundaries can be obtained by downloading publicly available data on tidal harmonic constants covering most of the oceans. TPXO (Egbert et al., 1994; Egbert and Erofeeva, 2002) is a widely used global tidal data. It is a series of fully global models of ocean tides, which best fits, in a least squares sense, the Laplace Tidal Equations and altimetry data. The TPXO models include eight primary (M2, S2, N2, K2, K1, O1, P1, Q1), two long period (Mf, Mm) and 3 non-linear (M4, MS4, MN4) harmonic constituents (plus 2N2 and S1 for TPXO9 only). More detailed information can be found at 🔗2). In addition, also the NAO.99b tide model (Matsumoto et al., 2000 and 2001), which is developed by assimilating TOPEX/POSEIDON Altimeter Data into Hydrodynamical Model, can provide global tide data. This model provides 16 short-period harmonic constituents (M2, S2, K1, O1, N1, N2, P1, K2, Q1, M1, J1, OO1, 2N2, Mu2, Nu2, L2, L2, and T2), 7 long-period harmonic constituents (Mf, Mm, Msf, Msm, Mtm, Sa, Ssa) data with of 0.5 degrees, and provides 16 short-period harmonic constants around Japan with a resolution of 5 minutes. More detailed information can be found at 🔗3.

Once the model and tidal boundaries have been established, they need to be tuned and calibrated before using. The model can be run with tidal boundaries and tidal potential. The tidal results are more sensitive to changes in the bottom friction coefficient and the bathymetry. By calibrating the tidal simulation, a reasonable bottom friction coefficient can be set for the model. At the same time, the difference in the tidal results caused by the bathymetry from dissimilar sources helps to find more suitable bathymetry data for the model.

In storm surge or storm surge flooding modelling, the input from meteorological forcing mainly includes surface wind shear stress and atmospheric pressure at sea surface level. In deep water, the sea level rises are mainly caused by the atmospheric pressure gradient, i.e. the water level rises approximately 1 cm at every reduction in pressure of 1 mbar. In shallow water, estuary and nearshore, wind shear stress is the dominant force of the storm surge, and the sea level rise is proportional to the square of wind speed. Therefore, accurate meteorological inputs, especially sea surface wind, is essential for storm surge modelling. The accuracy of storm surge results depends largely on the quality of meteorological data.

Depending on the storm surge modelling purposes, the types of meteorological forcing input are different. At present, there are two main sources of wind field for storm surge modelling: atmospheric model and empirical formula. Atmospheric models can provide global or regional meteorological fields; the main elements required for storm surge calculation are sea level pressure and 10 metres wind. The horizontal resolution of these data is between ten and tens of kilometres, and the forecast period can reach up to 240 hours. This kind of meteorological field is mainly used in the calculation of storm surge caused by extratropical cyclones. Compared with extratropical cyclones, tropical cyclones are smaller in scale but stronger in intensity, and atmospheric models cannot resolve the structure of tropical cyclones well. Therefore, the meteorological field from atmospheric models is not applicable to the typhoon storm surge calculation. Empirical formulas for tropical cyclone pressure and wind are often applied to create meteorological forcing fields for tropical storm surge models. Since the wind speed and pressure structure of tropical cyclones are close to axisymmetric, their distribution can be reasonably represented with an empirical formula for the radial distribution of wind or pressure.

The commonly used empirical formulas for pressure distribution mainly include the following:

Takahashi (1939):

Myers (1954):

Myers (1954):

Jelesnianski (1965):

Bjerknes (1921):

Holland (1980):

Where: 

• P_∞ = the ambient pressure;
• P_o = the pressure at the tropical cyclone centre;
• P(r) = the pressure at _r from tropical cyclone centre;
• A and B = empirical parameters that control the tropical cyclone size.

The tropical cyclone wind field is formed by the superposition of two vector fields (Ueno, 1981). The first vector field is a wind field symmetrical to the centre of the cyclone. The wind vector passes through the isobar and points to the left with a 20° deflection angle.

The wind speed is proportional to the gradient wind, which can be expressed by the following formula:

• V_r = the maximum wind speed;
• R = the radius of the maximum wind.

The second vector field caused by the movement of cyclone is superimposed on the wind system for the stationary symmetric cyclone, and that the wind velocity V→_f is:

• V→_t = velocity of typhoon;
• R = the radius of the maximum wind.

Consequently, the wind velocity W→ is:

The typhoon centre pressure, maximum wind speed, moving speed and other parameters can be obtained from websites of national meteorological agencies.

The governing equations of numerical simulation were determined by the flow dynamic theory.

The three-dimensional flow equations are as follows:

• Continuity equation:

• Momentum equation:

Based on hydrostatic approximations and incompressible assumption of fluid, the depth integrated two-dimensional storm surge governing equations can be written as:

• Continuity equation:

• Momentum equation:

Where:

• ξ= free surface elevation relative to the geoid;
• h = water depth;
• f = Coriolis coefficient;
• ρ = density of water;
• g = gravitational acceleration;
• (τ_x,S , τ_x,b) = free-surface shear stress in x and y direction;
  
• W_ = wind speed at 10 metres above sea surface;
• W_x , W_y = wind speed components in x and y direction;
• C_d = wind Drag coefficient which is relevant to wind speed;
• _{τx,b, τy,b} = bottom shear stress in x and y direction;
  
• u, v = depth-averaged horizontal velocity components in x and y direction;
  
• P = atmospheric pressure at the free surface;
• ε = depth-averaged horizontal eddy viscosity.

Hydrodynamic models are generally divided into vertically averaged 2D barotropic, 3D barotropic, or 3D baroclinic models. These types of models can be used for storm surge modelling. Minato (1998) first studied the effect of a 3D model on storm surge results and simulated the water level change caused by Typhoon 7010 in Tosa Bay, Japan. The results showed that the difference between the 3D model and the 2D model is about 2 to 10%. Weisberg and Zheng (2008) found that under the condition of setting the same bottom friction coefficient, the 3D model simulates higher extreme water level than the 2D model. Subsequently, simulating the storm surge of Hurricane Ike in the Gulf of Mexico, Zheng et al. (2013) pointed out that the 2D and the 3D models have differences in the trend and peak values of water level, but the calibration of the bottom friction is more important. Ye et al. (2020) studied the influence of baroclinic models on storm surge simulation results. Sensitivity tests show that the impact of the baroclinic model on storm surge is not significant, but it has a greater impact on current.

Since the vertical velocity distribution structure of the 2D model is different from that of the 3D model, the vertical average velocity is greater than that near the bottom layer, so that the bottom shear stress of the 2D model will be greater. Satisfactory results can be obtained for both types of models by calibrating the bottom friction.

Regarding operational storm surge modelling, computational efficiency is a factor that must be considered. 3D models generally divide the water body into multiple layers vertically, so they require more computing time than 2D models. If the purpose of storm surge modelling is to obtain water level, rather than currents, the 2D models are the best choice for balancing calculation efficiency and accuracy. Most of the existing operational storm surge forecasting systems are 2D models.

During the simulation of storm surge inundation, the grid points close to the coastline in the model will be either wet or dry due to the fluctuation of water level. Therefore, the model needs to determine the dry and wet state of a grid point according to the state of the surrounding grid points.

Assuming that the state of a certain grid point is wet, if the calculated water level makes the water volume less than zero, then this grid point will become a dry grid point and will not participate in the next calculation. In practical calculations, to prevent the negative value of the water volume at a grid point, which makes the momentum equation and continuity equation meaningless, a threshold value greater than zero is usually selected. When the water volume is less than the threshold value, the state of the grid point is defined as dry.

Assuming that the state of a certain grid point is dry, the first step is to check how many of the surrounding grid points are wet grid points. If more than one is wet grid point, then the water level is averaged over these wet grid points. If the averaged water level is greater than the threshold water level, the state of this grid point may become wet, otherwise it still remains as a dry one. In the second step, it is necessary to further check the transport over the cross-section area between the grid point and the surrounding wet grid points. If these transport cross-section areas are all positive, then the dry grid point becomes a wet grid point and participates in the next calculation; otherwise, it is still a dry grid point.

The grids used in most sea level models are mainly divided into two categories: structured and unstructured grids. The structured grid nodes are arranged in an orderly manner, and the connection relationship between adjacent nodes is fixed. In contrast, the unstructured grid nodes are arranged in an unordered manner and the adjacent nodes have no fixed connection relationship. Differently from structured grids, unstructured grids with triangular elements allow to adjust the resolution flexibly to depict complex shapes of coastline and estuary (Figure 7.10).

Structured grids mainly have two forms: rectangular grids and curvilinear grids. Relatively speaking, the curvilinear grid can adapt better to the complex shape of the coastline, and it can also realise the change of grid resolution that is more advantageous for storm surge simulation in estuary area (Figure 7.11).

Triangular grids were the main forms of unstructured grids until a new form of unstructured grids, the SCVTs (Figure 7.12), appeared in ocean modelling a decade ago (Ringler et al., 2013). Like the triangular grids, the SCVTs can adapt well to complex coastline, and facilitate a smooth transition from coarse resolution grid cell to high resolution grid cell. In addition, they also solve the computational instability problem caused by small acute angles in the triangular grid. At present, the SCVTs model has been used in the China Sea (as shown in Figure 7.13)

The main discretization methods include finite difference method, finite element method, and finite volume method. Discretization methods also correspond to the types of grids used. The finite difference method is generally used for structured grids, while the finite element method and finite volume method are generally used for unstructured grids.

The FDM is one of the simplest and oldest methods to solve the storm surge problems, and it is still widely used. This method divides the solution domain into differential grids, replacing the continuous solution domain with a finite number of grid nodes. By using the Taylor series expansion, the derivative of the governing equation is discretized by the difference quotient of the function value on the grid node, so as to establish the algebraic equations with the value on the grid node as the unknown quantity. This method is an approximate numerical solution that directly transforms a differential problem into an algebraic problem. The mathematical concept is intuitive and simple to express. It is an earlier and relatively mature numerical method.

The basic difference expression mainly has three forms (Figure 7.14): forward difference, backward difference, and centre difference. The first two formats are first-order derivatives, while the last format is second-order derivatives. Different computational schemes can be obtained through the combination of several different formats of time and space.

According to the precision of the scheme, we can distinguish among first-order, second-order, and high-order accuracy schemes. Depending on the influence of the time factor, the difference scheme can also be divided into explicit scheme, implicit scheme, explicit and implicit alternate scheme.

The FEM is based on variational principle and weighted residual method. The basic solution idea is to divide the computational domain into a finite number of non-overlapping elements. In each element, some appropriate nodes are selected as interpolation points of the solution function. The variable in the differential equation is changed into a linear expression composed of the node value of each variable or its derivative and the selected interpolation function. The differential equation is solved discretely by means of variational principle or weighted residual method.

According to the difference of the weight function and the interpolation function, the finite element method is divided into several computational schemes. For the choice of weight function, there are collocation methods, moment method, least square method, and Galerkin method. According to the shape of the computing cell grid, there are triangular grid, quadrilateral grid, and polygonal grid. Triangular grids are commonly used in storm surge modelling. The accuracy of the interpolation function is divided into linear interpolation functions and high-order interpolation functions. Different combinations also constitute different finite element calculation schemes.

The FVM is also called the control volume method. The basic idea is to divide the computational domain into a series of non-overlapping control volumes, and make a control volume around each grid point. Integrate the differential equations to be solved for each control volume to obtain a set of discrete equations. The unknown is the value of the dependent variable at the grid point. In order to obtain the integral of the control volume, it is necessary to assume the changing law of the value between grid points, i.e. the distribution profile (continuous or segmented) of the assumed value. From the selection method of the integral region, the finite volume method belongs to the subregion method in the weighted residual method. From the approximate method of the unknown solution, the finite volume method is a discrete method using local approximation. The physical meaning of the discrete equation is the conservation principle of the dependent variable in a finite controlled volume, just as the differential equation expresses the conservation principle of the dependent variable in an infinitely small controlled volume. The discrete equation obtained by the finite volume method requires the integral conservation of the dependent variable to be satisfied for any set of control volumes, and naturally also for the entire computation area. This is the attractive advantage of the finite volume method. Some discrete methods, such as the finite difference method, only satisfy the integral conservation if the grid is extremely fine. The finite volume method shows accurate integral conservation even in the case of coarse grids. As far as the discrete method is concerned, the finite volume method can be regarded as an intermediate between the finite element method and the finite difference method.

Numerical simulation of storm surges began in the 1950s. After decades of development, it has emerged that many numerical models can be used to storm surge simulations. Commercial models include MIKE21 and TuFlow, while examples of free models are: ADCIRC, Delft3D-FLOW, POM, FVCOM, ROMS, and SCHISM. Free numerical models generally provide the source code of the model, so that the model can be modified as needed when establishing a forecasting system. The models listed in Table 7.1.can be used to establish a complete operational storm surge forecasting system.

Table 7.1. Geometric interpretation of difference expression.

Data assimilation techniques are used to combine model and observed data to obtain the best estimate of the state of a system (see Chapter 4.3 for more details). Statistical techniques are often employed to find a solution which, ideally, minimises some error metric. For storm surges, this is done to obtain fields of sea surface height that can help us to better understand past events or to improve the quality of forecasts. An overview of the application of data assimilation to storm surge modelling and forecasting is provided in this section. Henceforth, references to errors mean some metric distance between the model and observations.

In order to reduce errors in storm surge models, especially for forecasting in which the accuracy may influence real-time decision making, it is important to understand the sources of error. Some of the main sources are given below.

Quality of input datasets

This includes atmospheric surface forcing and tidal forcing at the boundaries (see Section 7.2.2). Storm surges are largely forced phenomena; therefore, the accuracy of forcing is key and errors in the related datasets may be transferred into the storm surge component of the modelled sea surface height. Errors in input datasets may arise from similar sources as the storm surge model, including model and instrument errors.

For example, errors in tidal amplitudes and phases at the boundaries will propagate with the tidal waves into the domain. As discussed in Section 7.2.1, interactions with the tides can influence both timing and height of a storm surge, so it is important to have an accurate tidal component in the model. The accuracy of the surface atmospheric forcing is important as well, especially the components of wind and surface pressure. Due to the forced nature of storm surges, this is one of the largest sources of error in storm surge forecasts (Horsburgh et al., 2011).

Tuning of model parameters

It is common practice to adjust various model parameters to obtain a better solution. For example, it could include the tuning of bottom/surface friction coefficients. It is very unlikely to find a perfect parameter set, and the iterative processes often used can lead to non-optimal solutions. See Section 7.2 for more information about parameters used in storm surge models.

Representativity errors

Representativity errors arise from models’ ability to represent variables and processes such as resolution (Daley, 1991). For example, a coarse model may not be able to resolve finer scale features, which is present in the observations. For storm surges, this is particularly important nearby complex coastlines and estuaries. Similarly, coarser models may mean smoother bathymetry in these areas, which can significantly affect the modelled surge.

A model may not simulate all processes required to accurately model a storm surge. For example, if tides are not included in the model, only the atmospherically forced component of sea surface height is being generated, and contributions from tide-surge interactions will be missing. Other examples include the lack of tidal processes such as self-attraction and loading, or not including the inverse barometer effect in the model.

An attempt to reduce the impact of the errors described in the previous section can be made using data assimilation. For storm surge modelling, assimilation of observations may occur directly into the model or indirectly via input datasets.

Datasets used as atmospheric forcing often contain assimilated observations. The generation of the storm surge is highly dependent on the model’s interaction with these datasets and it is vital that they are accurate. For example, the forecasted atmospheric fields used at the UK Met Office use assimilation of atmospheric observations. There are also many reanalysis datasets available, such as the ECMWF ERA5 dataset (Hersbach et al., 2020) that assimilates observations from multiple sources to generate atmospheric data. While these examples are suitable for extratropical storm surges, they may not sufficiently resolve intense tropical cyclones, meaning that parametric methods may be a better option (see Section 7.2.2.3). There are also assimilative alternatives, such as the MTCSWA datasets (Knaff et al., 2011), which blend together parametric fields and observations. These have been shown to have some benefit for forecasting (Byrne et al., 2017). The same is true for the datasets used to derive tidal signals at the model boundaries. Examples of such datasets include TPXO (Egbert et al., 2002) and FES (Lyard et al., 2021), which incorporate data from satellite altimetry and tide gauges. See Section 7.2.2.2 for more information on these tidal datasets.

Sea surface height may also be assimilated directly into the modelled sea surface. There are two sources used, both with advantages and disadvantages: tide gauges and satellite altimetry. Tide gauges (and other fixed instruments such as bottom pressure recorders) offer information that is frequent and consistent in time, making them useful for capturing ocean processes of all frequencies (including storm surges). However, they are generally spatially sparse. On the other hand, altimetry data offer good spatial coverage but are less consistent in time, as a satellite only returning to the same location once every number of days. This makes the data useful for longer periods of periodic ocean processes.

Tide gauge data are currently assimilated for storm surge forecasting at some institutions (see Section 7.2.4.4). There are representativity challenges that must be considered when using these data. Most importantly, modelled sea surface height variables and observed variables must represent the same physical quantity. For example, do both datasets contain the same components of sea level such as tides and inverse barometer effect? The datum on which the data are based must also be considered. The sea level anomaly can be used to overcome these problems if a long enough record is available.

For real time forecasting, data assimilation is used to generate an improved initial condition for a forecast model run. The forward propagation of errors can be reduced by creating a more realistic initial condition. This is important to improve the lead times over which good forecasts may be given. The use of data assimilation for storm surge forecasting has been shown to offer improvements over short lead times (Heemink, 1986; Verlaan et al., 2005; Madsen al., 2015; Zijl et al., 2015; Byrne, 2021). However, the duration of these improvements may be limited to a few hours of forecast, due to the forced nature of storm surges.

This improved initial condition can be generated by running the model for some historical period up until today, including forcing with assimilated observations and potentially direct assimilation of sea surface height. This model run can effectively be seen as a continuous simulation with assimilative steps at some predefined frequency, for example every 6 hours. When a forecast is desired, the most recent state from this model can be taken and used as the initial condition. A forecast simulation is then done using no assimilation, as no observations are available in the future. This means that the atmospheric forcing used is also a forecast. Figure 7.15 illustrates this process.

There are several methods that have been used with success in storm surge modelling, including Optimal Interpolation (Gandin, 1966; Lorenc, 1981 and 1986; Daley, 1991), variational assimilation (Lorenc,1986), Kalman filters (Kalman,1960) and Ensemble Kalman Filters, (Evensen, 2004). In all methods, a key step is the estimation of spatial error covariances in both the model and the observations. This can be parametrically, as shown in the example in Figure 7.16, or by deriving covariances from an ensemble of model states. An example of the latter is described in Section 7.2.4.4.

Data assimilation has the potential to add significantly to the computation and time resources required, especially for ensemble systems. As noted, for real-time forecasting systems is vital that a balance is made between accuracy and speed, i.e. useful forecasts need to be delivered in a timely manner (Horsburgh et al., 2011).

The examples below are correct at the time of writing.

UKMO

The UK Met Office provides storm surge forecasting for the United Kingdom. Its 2D operational model does not currently assimilate any data into the model sea surface height. The atmospheric forcing used does include assimilated observations and comes from UKMO or the European Centre for Medium-Range Weather Forecasts (ECMWF) models, depending on the system. They also have more general operational 3D models that assimilate temperature, salinity and sea level anomaly, but the last is only done in deep water.

Rijkswaterstaat

Rijkswaterstaat provides storm surge forecasts for The Netherlands. Its system assimilates information from tide gauges around the Northwest European Shelf, especially in the North Sea (Verlaan et al., 2005; Zijl et al., 2015). They use a steadystate Kalman Filter (SSKF), which uses a stationary Kalman gain derived from an ensemble of states, such as might be used for the ensemble Kalman filter (EnKF). SSKF offers more computational efficiency than EnKF, and potentially a better representation of the error covariance than the standard Kalman Filter. Additional localization steps are also applied to the assimilation, to limit the distance from observations over which information is assimilated.

Like any other forecasts, sea level predictions have an associated uncertainty. The threat to life and property of extreme sea level events makes estimation of this uncertainty and the generation of a range of possible water levels (probabilistic forecast) particularly important for risk managers and decision-makers. The error of a single forecast time series can be assessed by comparison with in-situ tide gauges at specific locations and grid points. However, uncertainty of the forecast and its dependence on the forcing, model characteristics, and set up is usually unknown.

Uncertainty of sea level forecasts depends on several factors and may contain errors in both the tide and the non-tidal residuals (storm surge) components. During a storm, the storm surge is mainly driven by the weather conditions at the sea surface. This is considered to be the dominant source of uncertainty in sea level forecasts and may change significantly depending on the meteorological conditions. For this reason, ensemble storm surge forecasts based on weather ensemble prediction systems (EPS) are the most common approach to generate probabilistic forecasts (forecast plus a confidence interval).

The weather is a chaotic system highly sensitive to the initial state (Lorenz, 1965) that can only be deterministically predicted for about 10 days. Therefore, the standard procedure for dealing with forecasts uncertainty, i.e. the combination of different model solutions or ensemble modelling, was initially applied to meteorological forecasts (Leith, 1974; Hamill et al., 2000). Conceptual background of ensembles is chaos theory; they are a valuable tool to deal with equations in which several nonlinear processes and interacting variables are present. This is the case of meteorological models but also of ocean models and, particularly, of tide and surge models. Hence, their application is today strongly recommended in oceanography.

An EPS is based on the combination of a set of forecasts with different controlled changes in the initial conditions, the model physics or the open boundary conditions (Palmer and Williams, 2010; Flowerdew et al., 2010). All these modifications are designed to represent the uncertainties in the knowledge of the weather state. For example, different initial conditions allow to include those perturbations that grow most rapidly in time, in a context where slight changes to the initial conditions may lead to significantly different forecasts (Buizza and Palmer, 1995). Slight modifications of the set of equations (including different values in the parameterization constants representing different processes) also provide estimation of model uncertainties contributing to the forecast error.

Deviation of wind and sea level pressure fields from their actual evolution will determine a corresponding deviation on the predicted sea level. Therefore, their uncertainty derived from weather EPS will cause an uncertainty in the evolution of the sea level, linked to the meteorological forcing and affecting mainly the storm surge component. If different weather predictions are used to drive different sea level simulations, the probability distribution function of the forecast sea level values allows estimating the uncertainty of the sea level forecast and the probability of exceeding a given sea level threshold.

The ECMWF EPS has been operational since 1992 (Molteni et al., 1996). It was first applied to storm surge operational forecasts in the North Sea by Flowerdew et al. (2010 and 2013), who provided skilled probabilistic forecasts of sea level and showed that ensemble spread was a reliable indicator of uncertainty during large surge events.

Storm surge ensemble predictions have been used to forecast sea level in Venice by Mel and Lionello (2014a). They used a 50 members ensemble to simulate 10 events showing that EPS slightly increases the accuracy of the prediction with respect to the deterministic forecast, and that the probability distribution of maximum sea level produced by the EPS is acceptably realistic. They also showed that the storm surge peaks correspond to maxima of uncertainty and that the uncertainty of such maxima increases linearly with the forecast range. The same procedure was used by Mel and Lionello (2014b) for the simulation of the operational forecast practice for a three-month period (fall 2010). It revealed that uncertainty for short and long lead times of the forecast is mainly caused by the uncertainty of the initial condition and of the meteorological forcing, respectively. The probability forecast based on this ensemble technique has a clear skill in predicting the actual probability distribution of sea level. A computationally cheap alternative, called ensemble dressing method, has been proposed by Mel and Lionello (2016). It replaces the explicit computation of uncertainty by ensemble forecast with an empirical estimate. Instead of performing multiple forecasts, the procedure “dresses” the forecast of sea level with an error distribution form, which includes, on one hand, a dependence of the uncertainty on surge level and lead time and, on the other hand, of the uncertainty of the meteorological forcing. This computationally cheap alternative also provides acceptably realistic results.

Apart from the meteorological input, other sources of error on sea level forecasts can be attributed to the ocean model characteristics and/or to the setup of the system: bathymetry, spatial resolution, model domain, tidal forcing, temporal resolution of the meteorological input, barotropic or baroclinic models, ocean open boundary conditions, etc. Currently, sea level variations on timescales of hours/days are operationally forecasted through different barotropic and baroclinic models, sometimes over the same area. Therefore, another option is the combination of existing operational models with different characteristics, forcings and even physics (multi-model forecast).

A multi-model storm surge forecast was first implemented by Deltares (an independent Dutch institute) in 2008, combining existing operational storm surge forecasts from different countries in the North Sea. The system included the use of the Bayesian Model Average (BMA) statistical technique for validation of the different members and generation of a combined improved prediction, with a confidence interval (Beckers et al., 2008). In the same year, this methodology was tested for the Spanish coast by Puertos del Estado (Spain) (Pérez et al., 2012), combining the output of Nivmar (Álvarez-Fanjul et al., 2001), an existing storm surge forecasting system, with circulation (baroclinic) models already operating in the region. Nowadays, at Puertos del Estado is operational a multi-model surge forecast named ENSURF that combines Nivmar with two Copernicus Marine Service regional operational models, IBI-MFC (Sotillo et al., 2015) and MedFS (Clementi et al., 2021).

The BMA technique requires near-real time access to tide gauge data and automatic quality control of this data (as required by the Nivmar system as well), and specific data tailoring of model outputs. It is applied to the surge or non-tidal residual component of sea level because this can be approximated by a normal distribution (which is not the case for total sea level including tides, especially for strong semidiurnal regimes). So, observations from tide gauges and model data for those models providing total sea level must be previously decided. This could be considered a limitation but, in practice, it is the best way of optimising the final total sea level forecast by using the tidal component obtained from historical tide gauge observations at each site. ENSURF is also a valuable operational validation tool that allows a detailed assessment of the skills of different models to forecast coastal sea levels. A first deterministic forecast is provided by the old Nivmar solution early in the morning every day and, when later the Copernicus Marine Service forecasts become available, they are integrated with the tide gauge data and, by means of the BMA technique, a probabilistic forecast band is generated for each harbour (Pérez-González et al., 2017, Pérez-Gómez et al., 2019) (Figure 7.17).

Generally, use of the probabilistic methodology improves the forecast and gives significant added value to existing operational systems, as there is no single model that outperforms at all tidal stations and synoptic conditions. However, further work must be done with the BMA technique to predict the storm peaks which, in some weather conditions, are better captured by single systems.

A multi-model ensemble forecasting system has been recently developed for the Adriatic Sea combining 10 models predicting sea level height (either storm surge or total water level) and 9 predicting waves characteristics (Ferrarin et al., 2020). Other examples of this technique can be found in New York (Di Liberto et al., 2011) and the North Sea (Siek and Solomatine, 2011).

Storm surge models have been traditionally validated with time series of coastal sea level measured by tide gauges. These data allow assessing the skills of the model to reproduce observed water heights at specific points along the coast. Note the advantages this application presents with respect to sea level data from satellite altimetry, less reliable along the coastal strip and with a lower temporal resolution. Fortunately, there are hundreds of tide gauges around the world that become a very valuable and reliable dataset for storm surge validation (Muis et al., 2016 and 2020, Fernández-Montblanc et al., 2020). In some cases, these stations provide ancillary meteorological data, such as wind and atmospheric pressure, which can also be used to validate the model meteorological forcing. In addition, tide gauges transmit data in near-real time that can be integrated in an operational validation of the forecasting system (Álvarez-Fanjul et al., 2001).

In most cases, the forecast will provide the tide and storm surge signals (hourly to daily timescales), usually dominant at the tide gauge records, but will not be able to reproduce higher-frequency sea level oscillations such as infragravity waves, seiches or meteotsunamis, with periods of the order of a few minutes. It is important to define which observed “sea level” product will best fit the validation purpose, according to the physical processes included in the system. The most adequate standard product for existing operational storm surge forecasting systems are filtered hourly values from tide gauges. As new models include additional processes (e.g. fully coupled models including wave effects), higher resolution bathymetries, and forcing fields, the use of lower temporal sampling data will become more important and the validation process more challenging.

Normally, the model output at the grid point closer to the tide gauge is selected. However, the validation results will depend on the resolution and quality of the bathymetry data, as well as on the location of the tide gauge: if it is in an open site or inside a harbour or bay with important local effects, it may not be resolved by the model.

A careful validation of both tide and surge components should be performed to verify not only the final total sea level forecast, but also the quality of the tidal signal in the model, which can be an important source of error especially on shallow waters with high tidal range. For these reasons, model and tide gauge data must be de-tided, applying a harmonic analysis to both time series, as well as computing the tide prediction for the analysed period, and the surge or non-tidal residual after tide subtraction. The performance assessment can then be made in terms of comparison of the harmonic constants (amplitude and phase) from model and observations, and in terms of model data time series comparison of tide, surge and total sea level.

It is important to mention that sea levels measured by the tide gauge will be related to a local, regional or national datum. The model forecast is theoretically referred to mean sea level, though this mean sea level may be affected by the model setup implementation, domain, etc. Therefore, the mean should be subtracted from both time series at each location before comparison.

Metrics for time series validation can be found at Section 4.5.1. Most of these metrics describe the overall performance of the model for a time period of several days (the storm duration), months or years. However, they do not reflect the predictive skill for extreme surge events. This can be better evaluated, for example, in terms of the differences in the highest percentiles (e.g.: 95th, 99th percentiles) or the maximum observed and modelled value. For validation of long time series (multi-decadal hindcasts), it is possible to use annual maxima, annual percentiles, and extreme sea levels for specific return periods, obtained through extreme sea level analysis (Muis et al., 2020).

Taylor diagrams can be used to graphically indicate the performance of different competing models or solutions, providing information of the Pearson correlation coefficient, the standard deviation and the RMSE at each tide gauge.

Progressively, the storm surge models will consider inundation, and additional validation of the flooding extent will be required. This is less straightforward and requires other types of data, such as locations of flooded points, marks left by the water or reports about the flood chronology (Le Roy et al., 2015). This information is commonly available after the event for a delayed mode validation; e.g. for validation of inundation, Loftis et al. (2017) used crowdsourced GPS data and maps of flooded areas obtained by drones.

The main outputs of storm surge models are: time series output, maximum elevation field output (extreme values at every time step for water surface elevation), ensemble forecast elevation field output, animation output.

The time series output is usually plotted in a two-dimensional rectangular coordinate system, the abscissa is time and the ordinate is water level. The time series output of storm surge models are the water level changes at a certain location. In order to facilitate the comparison between the results of simulations and the observation data, multiple result curves can be plotted in the same coordinate system.

Generally, the results of storm surge models (without tide) can be directly used for plotting time series diagrams, but sometimes attention should be paid to the change of the total water level at a certain point. Therefore, there is a way to output the total water level, that is the results of the storm surge model without astronomical tide directly superimposed on the astronomical tide from harmonic analysis, obtaining in this way the time series results of the total water level. Figures 7.18 and 7.19 show examples of time series model storm surge result and storm surge superimpose on harmonic analysis tide.

Among the output methods of storm surge models, there is an output form called maximum elevation field. This kind of field output is not the water level field at a certain time, but extracts the highest water level value of each grid point in the calculation process to form the maximum elevation field. The maximum water level field can be used to grasp the distribution of the maximum water level during a Typhoon and to identify the more severely affected areas along the coast. Figure 7.20 shows the maximum storm surge field during 2019 Typhoon Mitag (1918).

The storm surge ensemble forecast usually uses the respective meteorological forcing fields of the ensemble members to calculate the storm surge separately. The output of ensemble forecast fields mainly includes the following forms: (a) ensemble mean field; (b) probability field; and (c) postage stamp maps.

a. Ensemble mean storm surge field

In order to obtain a definite forecast result, it is necessary to synthesise the respective results of the ensemble members. It is generally used to assign different weights to the results of each member and to superimpose the results of all members, i.e. the weighted average method (Wang et al., 2010). The superimposed result is output in the form of the elevation field, and the ensemble forecast field is obtained as result.

The track map of Typhoon Mitag can be seen in Figure 7.21. The storm surge results of the subjective typhoon forecast track, fast track, slow track, left track, and right track are used to synthesise the ensemble forecast water level field applying the weighted average method (Figure 7.22). In this example, the weight of the storm surge result of the subjective forecast typhoon track is 60%, while the weight of the storm surge result of the other tracks are all 10%.

b. Probability storm surge field

The typhoon ensemble forecasting tracks for storm surge numerical simulation can describe the surge field under different typhoon track scenarios. By equitably assigning weights to different track results, the probability field distributions under different extreme values of storm surges can be clearly displayed, and the intensity probability of coastal storm surges can be more intuitively presented (Liu et al., 2020). Figures from 7.23 to 7.25 show the probability distribution of forecast storm surge over 0.5m, 1.0m and 2.0m of Typhoon Mitag.

c. Postage stamp maps

A postage stamp map is a set of small storm surge maps drawn by the results of the individual members (Figure 7.26). Forecasters can learn about the possible situation of each ensemble member through the postage stamp map, thereby estimating the magnitude of the maximum storm surge and the range of impact. (WMO, 2012).

The output of the water level field is very helpful for grasping the distribution of the storm surge process over a whole region. The output of the storm surge water level field is the elevation value of all grid points at a certain time. The elevation field figure at each moment is taken as a frame of the animation, and all the frames are connected to form the elevation field animation (Figure 7.27). The elevation field animation can intuitively reflect the changes in the water level of the entire area during the impact of the storm surge.

After decades of development of storm surge numerical models, many countries have established their own operational storm surge forecasting models. For example, in the United States a storm surge forecasting system is operating through the SLOSH model, in which the wind field is established based on the cyclone path, maximum wind speed radius, storm centre, and environmental pressure difference; it provides operational forecast products and storm surge inundation guidance products (Jarvinen and Lawrence,1985). China established ver3.0 of the PMOST forecasting system, which is based on a depth-averaged two-dimensional shallow water equation in the vector invariant form, and uses a SCVTs grid. It can enhance the resolution in key areas and fit the coastline. The system is able to couple astronomical tides and simulate flooding processes. With the GPU acceleration technology, the efficiency of storm surge simulation along the coast of China can reach 60 sec/day. The system can also perform ensemble forecasts based on multiple typhoon events and storm surge probability forecasting. The Indian Institute of Technology storm surge model was developed in the 1980s and applied to storm surge forecasts in the Indian Ocean and the Arabian Sea. It uses rectangular Cartesian coordinates and separates land and water during calculations. It has been applied throughout the north Indian Ocean (Dube et al., 1984, Dube et al., 1985).

The two-dimensional storm surge model developed by the Japan Meteorological Agency also uses rectangular coordinates (Hasegawa et al., 2015). In the numerical calculation, the water and land are separated with the flexible mesh, i.e. fine grid is used in shallow water and coarse grid is used in deep water. The system can also provide ensemble forecasting products. In the mid-1980s, the Netherlands developed a numerical fluid dynamics model called the DCSM, which uses a depth integrated shallow water equation. The driving force of the model is provided by a high-resolution regional meteorological model. In the early 1990s, the Kalman filtering method was used in DCSM to assimilate the water level (Verlaan et al., 2005; de Vries, 2009). The UKMO developed the storm surge forecasting model CS3X, which is a tide-storm surge model. In this operational storm surge forecasting system, the interaction of tide and storm surge is considered. In recent years, a storm surge ensemble forecasting has been developed in this system (Flowerdew et al., 2013). A storm surge model covering the French overseas territories has been operated since the 1990s by Meteo-France (Daniel et al., 2009); it was established based on the spherical nonlinear shallow water equation. In order to solve the problem of the shore boundary, the C-grid difference format was adopted, with meteorological forcing provided by the Holland model (see Section 7.2.2.3). Table 7.2 provides a list and features of storm surge forecasting systems currently operating in various countries.

Table 7.2. List and features of operational storm surge forecast models.