Nivmar : A storm surge forecasting system for Spanish Waters *

The knowledge and prediction of sea level is important both for the scientific and economical point of view. To describe the circulation patterns in a region is convenient to characterize sea surface elevation. This variable is also a direct indicator of physical phenomena such a upwelling, very important for biological productivity. In the harbours daily work, the knowledge of sea level is critical (for example for dragging operations). When a storm surge is combined with severe waves, damage in the coastline structures can be important. All this reasons impulsed Puertos del Estado to develop a storm surge prediction system, called Nivmar (http://www.puertos.es/Nivmar), able to forecast and study the evolution of surface elevation. The evolution of sea surface elevation during a storm surge event is strongly influenced by bathymetry. In seas with a wide continental shelf, like the North Sea, surge evolution is dominated by the wind. In those seas, extra-tropical surges can be very large (up to 4 m in the North Sea). Along the Spanish coast, in the Atlantic Ocean (including the Canary Islands) and in the Mediterranean Sea (including the Balearic Islands), a common characteristic feature is found: the continental shelf is very narrow (Fig. 1). Therefore, residuals (difference between measured sea level and tidal elevation) are mostly produced by pressure variations. In fact, a first rough estimate of SCI. MAR., 65 (Suppl. 1): 145-154 SCIENTIA MARINA 2001


INTRODUCTION
The knowledge and prediction of sea level is important both for the scientific and economical point of view.To describe the circulation patterns in a region is convenient to characterize sea surface elevation.This variable is also a direct indicator of physical phenomena such a upwelling, very important for biological productivity.In the harbours daily work, the knowledge of sea level is critical (for example for dragging operations).When a storm surge is combined with severe waves, damage in the coastline structures can be important.All this reasons impulsed Puertos del Estado to develop a storm surge prediction system, called Nivmar (http://www.puertos.es/Nivmar),able to forecast and study the evolution of surface elevation.
The evolution of sea surface elevation during a storm surge event is strongly influenced by bathymetry.In seas with a wide continental shelf, like the North Sea, surge evolution is dominated by the wind.In those seas, extra-tropical surges can be very large (up to 4 m in the North Sea).Along the Spanish coast, in the Atlantic Ocean (including the Canary Islands) and in the Mediterranean Sea (including the Balearic Islands), a common characteristic feature is found: the continental shelf is very narrow (Fig. 1).Therefore, residuals (difference between measured sea level and tidal elevation) are mostly produced by pressure variations.In fact, a first rough estimate of SUMMARY: In this paper, a storm surge prediction system for the Spanish Waters is presented.The system, named Nivmar, is based on the ocean circulation Hamsom model and on the harmonical prediction of tides computed from data measured by the tide gauge network Redmar, managed by Puertos del Estado.Nivmar is executed twice a day, running Hamsom forced by meteorological fields derived from the INM (Instituto Nacional de Meteorología) operational application of Hirlam atmospheric model.Data from Redmar tide gauges is used to to forecast the tidal elevations, to validate the system and to perform data assimilation, correcting systematic errors in the mean sea level due to physicals processes that are not included in the ocean model (i.e. steric height).The forecast horizon is 48 hours.In order to validate the system with measured data from Redmar a very stormy 5 months period was selected.Results from this test (November 95 to March 96) are presented.Data from this experiment shown that Nivmar is able to correctly predict sea level in the region.A simple data assimilation scheme for sea level is described and results from its application are studied.Finally, special focus is made in future plans and potential developments and applications of the system.
the residuals in the Atlantic Spanish harbours can be obtained through the study of the time series of inverted atmospheric pressure records, although the effect of wind can be also of importance during severe storms.A special case is the Mediterranean Sea, where the inverse barometer effect can not be directly applied (Le Traon and Gauzelin, 1997) and the residuals are larger than the tides.Garret has developed a simple model of the Mediterranean (Candela et al., 1989;Le Traon and Gauzelin, 1997), which includes two basins and two straits (Gibraltar and Sicily).This model, driven by the mean atmospheric pressures on the two Mediterranean basins, produces better results than the simple inverse barometric correction and shows the important role of the Straits of Gibraltar when studying the barotropically induced sea surface elevation in the Mediterranean.
Results from the numerical simulations presented in this work confirm this critical role.In narrow shelf seas, the limited effect of the wind over the shelf produces smaller surges (up to 1 m).
All the operational surge prediction systems deals with the problem of correcting systematic errors in the mean sea level due to physicals processes that are not included in the ocean model (i.e. steric height).Nivmar is solving this problem by simply correcting the solution by the difference of the means in the previous year, but an assimilation scheme has been developed (section 4) and will be implemented in the future.The major problem (currently under solution) to set up the assimilation scheme is the reception in real time of the Redmar tide gauge network data.
This paper is structured as follows: first the Nivmar system is described making special focus on the ocean model implementation and on the technical details of the operational system.The next section is devoted to the validation of the system during a stormy period and, finally, the window based assimilation scheme, that will be implemented operationally in the future, is described and results from its application are shown.

THE NIVMAR STORM SURGE PREDICTION SYSTEM
The Nivmar System, consisting of a set of different applications and programs, will be used to predict sea level using the Hamsom model forced by atmospheric data generated by the INM (Instituto Nacional de Meteorología) with the Hirlam (HIgh Resolution Limited Area Model) atmospheric model (winds and pressures fields).The Hirlam forecasting system (Källen, 1996) is a limited area, primitive equation model developed by several European institutes and is based on hydrostatic approximation.
Hamsom is written in Cartesian coordinates.To account for the convergence of meridians, all distances in the horizontal axis are computed as a function of latitude and the cell volume so considered is distorted for mass conservation purposes.The code uses a two-time level numerical scheme (present and future).Some equation terms are treated semi-implicitly (like pressure gradient) and the remaining ones explicitly (momentum advection, horizontal diffusion).
Bottom stress is parameterised by a quadratic law in terms of current velocity: where Cb is the dimensionless drag coefficient.A classical value of 0.0025 was selected for all the runs.

Description of Nivmar
The system is based on the ocean circulation Hamsom model and on the harmonical prediction of tides computed from data measured by the tide gauge network Redmar.The forecast horizon is 48 hours.The model domain covers an area extending from 20ºN to 48ºN in latitude and from 34ºW to 30ºE in longitude (see Fig. 2).
When designing a storm surge prediction system it is important to have in mind that the barotropic set of equations used in the ocean model contain nonlinear terms.Those are the bottom friction (Eq.4), the advection of momentum (second and third term in Eqs. 1 and 2) and the product of elevation and current speed in the continuity equation (Eq.3).One of the effects of these terms is to transfer energy contained in one frequency to another (Bowden, 1983).For example, the M2 tidal constituent generates harmonics (M4, M6) with frequencies that are multiples of the original one.Also the residual elevation during a storm surge is modified by nonlinear interactions with the tide.In very shallow seas these nonlinear contributions cannot be neglected, becoming critical to predict the evolution of sea surface.In the North Sea and the English channel, for example, M4 becomes an important harmonic and residuals can't be properly predicted without the nonlinear contribution of energy from tides (Pugh, 1987).In those shelf seas, residuals are computed in the following way in order to take into account the nonlinear interactions between tidal and meteorological forcing: first, a complete simulation, including tidal and atmospheric forcing, is performed.Tides are introduced in the model by prescribing, in the open boundaries, the elevation induced by the harmonic constituents.For doing so, only a small subset of tidal components (typically 8) is usually available at the model boundaries.For computing nonlinear transfers of energy this is sufficient because those components contain most of the energy of the tide.Second, a tidal simulation of the same subset of harmonic components is carried out for the storm sim-ulation days.Finally, residuals are obtained by subtracting sea surface elevations from both simulations.Final elevation in selected coastal points is computed by adding those residuals to the tide prediction, which include all the astronomical components derived from harmonical analysis of measurements, not only the most energetic ones included in the simulations.
We have shown that, in narrow shelf seas (like the Spanish coast), where nonlinear transfers are not so critical, similar results to those computed by the previously explained method are obtained with the simpler method of computing residuals as the sea surface elevation derived from a simulation forced only by meteorological fields (Álvarez, 1998; Álvarez et al., 1998).Therefore, Nivmar is making forecasts of the residuals using this second method, which is simpler since only one model run is performed.The total sea surface elevation is computed, as in the first method, by adding the astronomical prediction, which includes all the constituents.In the case of Nivmar, this tide forecast is provided by the harmonical analysis of Redmar data.
Figure 3 is a scheme of the normal operations performed by Nivmar during one day.Nivmar is executed twice a day in cycles named 00H and 12H.The model needs input data (results from previous simulations for a warm start, meteorological fields from the INM application of Hirlam and sea level measurements) that are pre-processed at an initial stage.The model domain of the Hirlam model (resolution 0.5º x 0.5º) covers completely the Hamsom model domain and the data (6 hourly fields of winds at 10 m and pressures) are interpolated to the Hamsom grid by means of a bilinear interpolation routine (Press et al., 1992).The ocean model run covers a period spanning from 12 hours before the starting time (00 hours on the 00H cycle and 12 on the 12H cycle) to 48 hours later.Analysed meteorological fields (hindcasted fields computed using data assimilation) are employed on the first 12 hours of simulation.Results form Nivmar are placed on a massive storage system (Unitree) and distributed through the web.
Nivmar is under continuous development.The system started producing predictions in June 1998.Results from that year where used to create a first quality report (http://www.puertos.es/Nivmar/verif.html).During that period, the system showed a good behaviour in the Atlantic Ocean, but results in the Mediterranean Sea were not as good as expected due to a problem, already corrected, in the model domain (the eastern limit reached only 9ºE).In June 1999 the operational phase started for the Atlantic ocean.Results from 1999 with the extended bathymetry will be validated during the present year to make the system operational in the Mediterranean Sea.
The system is running at an X-Class Exemplar parallel computer at Puertos del Estado (Hamsom is running parallel in 4 HP PA-8000 processors).Updated predictions can be found under http://www.puertos.es/Nivmar/nivmareng.html.

Model implementation
The bathymetry employed in the simulations (Fig. 2), is based on the DTM5 data set (GETECH, 1995).It was built using a variable grid size scheme in order to reduce the number of computational points and, consequently, the computation time.This scheme is based on a zooming technique.A region spanning from 25ºN to 48ºN and from 20ºW to 30ºE keeps a constant resolution of 10' x 15'.The grid size in the rest of the domain is increased progressively towards the boundaries, being the size of each cell incremented by a factor of roughly 1.1 with respect to its inner neighbour.Small deviations are present in this factor in order to place the boundaries on a non-fractional latitude (20ºN) and longitude (34ºW).
In the North coast of the Iberian Peninsula, the available bathymetry did not include the existing narrow shelf.As a consequence, the effect of the wind was initially underestimated in the simulations.In order to improve the results of the model, the depth of the grid points near the coast were modified in this region.After some tuning, it was found that changing the depth of the two rows (or columns) of grid points adjacent to the land to depths of 10 m and 50 m produces optimal results.These changes in the bathymetry greatly improve the results of the model in the region.
Both in the benchmark presented in this paper and in the Nivmar system the time step employed is 10 minutes and the horizontal eddy viscosity is kept constant at means of the so called Charnock constant (α).After some tuning (Álvarez Fanjul et al., 1998), it was found convenient to use a high α value for the Charnock constant (α=0.032).The elevation at the Atlantic open boundaries is corrected by means of the inverse barometer effect.In the Mediterranean Sea the inverse barometer correction is not satisfactory (Le Traon and Gauzelin, 1997).Therefore, to model the whole basin would be desirable.Unfortu-nately, the atmospheric model only reaches 30ºE, leaving a small portion of the boundary conditions, including inverse barometer, the best solution was provided by placing a solid wall at the Eastern boundary.This is equivalent to reducing the area of the eastern Mediterranean basin and could have some impact on this region, but results presented in this paper show that good results are obtained at the western basin.

VALIDATION OF NIVMAR
The Nivmar system was validated with data from the Redmar tide gauge network (Pérez and Rodríguez, 1994).A very stormy period was chosen for validation (November 1995to March 1996).The largest surges ever recorded by Nivmar took place during this period.The atmospheric forcing employed in these test simulations (6 hourly fields of winds at 10 m and pressures) was obtained from the same source (fields supplied by the Spanish National Meteorological Institute produced by their operational system based on the Hirlam model) than the data operationally employed to force Nivmar.
Figure 4 shows a comparison of the residuals obtained by numerical simulation with those measured by the Redmar during the PROMISE Spanish Coast Data Set Period (see Fig. 1 for location of the tide gauges).The three stations of Figure 4 are representative of different regions of the Spanish Coast (Canary Islands, Mediterranean coast and Iberian Peninsula Atlantic coast).In order to establish a comparison, the frequency of the simulated data was changed from 10 minutes (time step of the model) to one hour (period between two measurements).The mean value of the simulated series was modified to the value of the measured one.These techniques are usually employed in hindcast studies of surge simulations (Vested et al., 1995).Table 1 shows a statistical comparison between simulated and measured series in the stations represented in Figure 1.
Model results in the Cantabric and Galician coast (stations from Bilbao to Vigo) are quite satisfactory, specially taking into account that this is a very stormy period.Maximum errors are around 20 cm, typical values for rmse are 5 cm and the correlation index is very high (0.95 for Vigo).Most of the largest peaks are correctly reproduced, including extreme storms (i.e.66 cm at La Coruña on January the sixth, corresponding to the 67 day of the simulated period).
Inspite of the presence of a narrow shelf, the effect of the wind on the sea surface elevation during severe storms cannot be neglected.A prediction computed by means of the Inverse Barometer Correction strongly underestimates the surge height dur-ing those events.Differences up to 43 cm between measured and computed residuals can be found, for example, in the large storm surge recorded at Bilbao on the 8th of February (the largest peak in Fig. 4).For this same event, the difference with the data generated from Nivmar was only 11 cm.This example shows the importance of taking wind into account in the simulations, even in a narrow shelf sea like the Cantabric.
The residuals measured at Bonanza show larger differences (rmax=69 cm) with the simulated data.This gauge station is located at the Mouth at the Guadalquivir river, and residuals are strongly influenced by the river outflow.Since this variable is not included in the simulations, the observed discrepancies are expectable.In periods where no important outflow from the river was present the agreement between measurements and simulations improves considerably.
At the Mediterranean sea the root mean square error and the maximum error have the same order of magnitude as in the Atlantic coast, having the correlation index lower values (minimum of 0.85 for Valencia).These lower values of the CI can be related to the quality of the winds.At the resolution of the atmosphere model (0.5º x 0.5º) wind solutions are better at the Atlantic coast than in the Mediterranean coast, because the eastward travelling atmospheric structures are affected by the complex topography of the Iberian Peninsula and, therefore, atmospheric fields are more complex to solve at the East side of the Peninsula (J.A. García-Moya, pers. comm., 1998).Future use of atmospheric fields of higher resolution (0.2º x 0.2º) will mean an improvement in the sea elevation solutions at this region.
It is interesting to note the important role of the Straits of Gibraltar in the evolution of sea surface NIVMAR: A STORM SURGE FORECASTING SYSTEM 151 elevation in the Mediterranean Sea. Figure 5 shows the results at Barcelona (compared again with measured residuals) of a similar simulation, but with the straits closed (the grid point corresponding to the straits was artificially converted to land).From the analysis of the figure it is clear that, if the Mediterranean sea is isolated from the Atlantic, the response of the model is worse.The statistical analysis confirms this idea: the correlation index goes down from 0.89 to 0.62, the rmse increases form 5.15 to 8.72 cm and the slope of the linear fit drops from 0.71 to 0.37.These results are coherent with the ideas behind the simple analytic model where subinertial flows through the straits and barotropically induced sea surface elevation in the two main Mediterranean basins are basically determined by the relations between the mean atmospheric pressures over the two basins and over the Atlantic (Candela et al., 1989;Le Traon and Gauzelin, 1997).
The agreement at the Canary Islands is highly satisfactory, where the rmse values are as low as 3 or 4 cm.

IMPLEMENTATION OF A VERY SIMPLE DATA ASSIMILATION
The window based data assimilation technique improves the results of the predictions by simply correcting the mean value of the simulated residuals.The correction is done by adding a constant value which is the difference of the mean values of the simulated and predicted time series during a lapse of time (window) in the immediate past (see Fig. 6).This simple technique, applied during the post-processing stage, helps to avoid systematic differences between the mean levels of the simulated and measured time series and greatly improves the quality of the solutions by reducing the bias.
The window must be large enough to filter the noise of the measurements (minimum of 1 day) and short enough to include data representative of the conditions existing in a near past instead of mean conditions (maximum 5 days).
This technique was applied to the simulated elevations generated in the 5 months validation period (see section 3) making use of the residuals measured by Redmar.Different combinations of window lengths (1, 3 and 5 days) and forecast horizons (6, 24 and 48 hours) were explored.Table 2 shows the results at the Bilbao station.As expected, the better results are associated with smaller forecast horizons (6 hours).Comparison with the results obtained before applying the assimilation shows the improvement obtained with this method.At Bilbao (for a forecast horizon of 24 hours and making use of a window of three days) better results are obtained for rmse (5.51 cm with data assimilation against 5.66 cm without assimilation), rmax (18.70 cm compared to 21.03 cm), m (0.89 to 0.82), b (0.79 to 1.44) and CI (0.92 to 0.90).Similar improvements were found in other locations.
It is also interesting to note that, for short forecast horizon predictions (6 hours), a small data window (1 day) is recommended, while for long forecast horizon (48 hours) a larger window (5 days) produces better results.

CONCLUSIONS AND FUTURE PLANS
A storm surge prediction system has been developed and validated for the Spanish coasts.Results of the system show a good behaviour and the dynamics of the Spanish coast, characterised by the presence of a narrow shelf (in contrast with other well studied regions like the North sea), are correctly reproduced.
Nivmar is a newly created system, but several applications and collaborations with other institutes derived from its data are already under way.At this moment, data from Nivmar are being used to study sea level evolution, together with altimeter data, both in the Mediterranean and the Atlantic, and the vertically integrated currents generated by the system were used in a project to monitor oil spills based on Synthetic Aperture Radar images (Envisys project).Results from the system are being also used, for example, by a shipyard (Naval Gijón) in order to forecast better time for shipping.
Several future plans, both for the short and long term, are already being developed: Short and medium term plans -Wind and tide induced surface current predictions.Hamsom will be used in tridimensional mode.
-Inclusion of the Nivmar results in the Puertos del Estado data base.
-Web based access to the Nivmar historical results (in graphical format).
-Application of a data assimilation scheme.
Long term plans -Introduction of temperature and salinity gradients in the simulations.This will allow Nivmar to predict and study those variables.
-Nested high resolution applications.Nivmar is already a useful system, both from the scientific and the economic point of view, but the potential of the system is much higher.TABLE 2. -Statistical comparison of measured and simulated residuals after applying the window data assimilation scheme, where win is the length of the data window in days, hori is the forecast horizon, Nr is the number of records,⎯ X is the mean value of the measured time series,⎯ Y is the mean value of the simulated time series after the assimilation, BIAS is the difference⎯ Y-⎯ X, rmse is the Root Mean Square Error, rmax is the maximum error, m and b are the slope and the interception of the linear fit and CI is the Correlation Index.
FIG.2.-Bathymetry employed in the numerical simulations.Resolution in the west part of the domain is poor due to the variable grid size scheme employed and to the fact that the first two model columns are equal in order to implement a zero gradient boundary condition for the transports.
FIG. 3. -Nivmar system box diagram showing the normal operations performed during one day of predictions.

TABLE 1 .
-Statistical comparison of measured and simulated residuals.Nr is the number of records, ⎯ X is the mean value of the measured residuals (mean value of simulated residual is corrected to this value), rmse is the Root Mean Square Error, rmax is the maximum error, m and b are the slope and the interception of the linear fit and CI is the Correlation Index.