Water circulation forecasting in Spanish harbours

1 Laboratori d’Enginyeria Marítima, Universitat Politècnica de Catalunya (LIM/UPC) and International Centre of Coastal Resources Research (CIIRC), Carrer Jordi Girona 1-3, 08034 Barcelona, Spain. E-mail: manel.grifoll@upc.edu 2 Institut Mediterrani d’Estudis Avançats (IMEDEA), UIB-CSIC, Carrer Miquel Marqués 21, 07190 Esporles, Spain. 3 Puertos del Estado, Av. Partenón 10, 28042 Madrid, Spain. 4 AZTI-Tecnalia, Marine Research Division, Herrera Kaia, Portualdea s/n, 20110 Pasaia, Spain.


INTRODUCTION
Anthropogenic pressure caused by economic harbour activities can alter the water quality inside harbour domains and affect adjacent areas that may be used for other activities such as fishing, swimming and nature conservation.Environmental protection policies aimed at fostering sustainable development have been applied across Europe.The EU Water Framework Directive (WFD, 2000/60/EC) has emerged as a common policy framework to protect continental, underground and coastal waters, and it also covers harbours and their adjacent waters.An understanding of water circulation patterns within harbours is an essential aspect in the management of issues related to water quality degradation (Hartnett and Nash 2004), risk analysis of pollution events (Grifoll et al. 2010) and accidental spills (Jordà et al. 2007).
Several operational products for harbour and coastal applications have been developed in recent years.Most of these applications are wave and water-level forecasting systems (see for instance, Carretero et al. 2000, Cox et al. 2002, Lin et al. 2008, or Fernández and Mayerle 2008).However, much less attention has been paid to water circulation products for harbour domains.Moreover, the experience gained in operational circulation systems implemented in the open sea (e.g.MFS system: Pinardi et al. 2003, or ESEOO system: Sotillo et al. 2007) may not be directly translated to harbour domains.Their reduced dimensions and intricate layout confer upon harbours restrictions which are not present in the open sea.In order to fill this gap, three Spanish harbours (Barcelona, Tarragona and Bilbao), Puertos del Estado (the Spanish National Ports and Harbours Authority) and the Maritime Engineering Laboratory of the Universitat Politècnica de Catalunya (LIM/UPC) signed a collaboration agreement.Among other objectives related to the monitoring of harbour waters through intensive field data campaigns, this agreement focuses on researching and modelling harbour hydrodynamic conditions with the aim of implementing an operational forecasting system able to provide daily forecasts of water circulation in harbours.Derivative products based on current forecasts, such as float trajectories, residence time maps and risk assessment of water quality degradation, are also of interest to this project.
At present, the system has been implemented in a pre-operational phase in the three aforementioned harbours (see locations in Fig. 1), and the results are only available for harbour managers.These three harbours are suitable test cases because of their different characteristics, and the learning gained from them may therefore be translated to most Spanish harbours.Barcelona and Tarragona harbours are located in the northwestern Mediterranean Sea, where tidal forcing is negligible, while Bilbao harbour is located in the southwestern corner of the Bay of Biscay, in a mesotidal environment.Barcelona harbour has two mouths (north and south, Fig. 1), while Tarragona and Bilbao harbours have only one.Finally, Tarragona and Bilbao have river discharges inside the harbour domain, while Barcelona does not.A first step in the implementation of any operational system is to evaluate the relative importance of the different forcings that are acting on the circulation, so that the system will be designed accordingly.This evaluation depends on the oceanographic characteristics of the harbours.Mestres et al. (2007) analyzed the Tarragona harbour hydrodynamics and found complex hydrodynamic patterns in the inner harbour domain, with a water circulation of a few cm s -1 .These authors suggest that wind stress is the main forcing of the water circulation in the harbour, inducing a predominantly barotropic circulation; they also emphasize that the baroclinic structure can have an important influence on the long-term circulation.The effects of the Francolí river, a small river with an averaged winter runoff of 1.4 m 3 s -1 that discharges into the harbour were also analyzed.The conclusion was that there were no significant effects of the river outflow on harbour circulation.In Barcelona harbour, the picture is different because there are two mouths.From an analysis of intensive field campaigns, Grifoll et al. (2011) found a predominant barotropic circulation, especially in winter.However, they concluded that wind forcing was not enough to explain all the current variability observed within the harbour.From numerical experiments, they showed that the existence of two mouths allows the shelf currents to play a significant role in the inner dynamics.Moreover, additional experiments demonstrated the minor influence of surface heat fluxes and freshwater sources.Bilbao harbour also shows a different picture due to the dominant influence of tides.Grifoll et al. (2009) carried out several numerical experiments to determine the relevance of different forcings on harbour circulation and found that tides were responsible for 60% to 90% of circulation variability, depending on the wind regime.They also found a moderate influence of the freshwater discharge from the Nervión river (a medium-sized river with an averaged winter runoff of 35 m 3 s -1 ) in the innermost part of the estuary.In summary, the environmental and morphological characteristics of each of the test harbours cause them to respond differently to the applied forcings.Therefore, the operational system configuration in each harbour may depend on these forcings.
This paper is devoted to presenting a Harbour Circulation Forecasting System (HCFS) developed for Spanish harbours.The aim is to describe the system and to justify the different choices adopted.An analysis of the sensitivity of the predictions quality to errors in the forcing fields is also addressed.Although the results presented here mainly focus on three harbours (Barcelona, Tarragona and Bilbao), the methodology and conclusions may be applied to most harbours, as we cover a wide range of possible configurations (microand mesotidal, one and two mouths, with and without freshwater discharges).The contribution is organized as follows: A brief description of the numerical codes and forcing fields is first given (System description).The experiments devoted to determining the optimal initialization strategy in each harbour (Initialization strategy) and a validation of the system results are presented next (System validation).The sensitivity of the forecast quality to forcing errors is then assessed (Sensitivity of harbour forecast to forcing errors).Following that, the results are discussed, underlining the main physical differences obtained among the harbours and the improvements required for the HCFS.Finally, the main conclusions are presented.

Nesting strategy
Inner harbour hydrodynamics may be influenced by processes occurring in the coastal vicinity.In consequence, a proper modelling of coastal processes outside the harbour domain is required in order to achieve good-quality results in the harbour.Moreover, coastal processes are in turn affected by the large-scale dynamics, so the modelling system must consider a wide range of scales, from the basin to the harbour scale.For this system we followed a strategy of using a hierarchy of one-way nested models, with four different levels of resolution (see Fig. 1).The first one is a regional model able to solve the mesoscale dynamics in the area (fronts, meanders, eddies, etc.).Thus, it covers a large area with a typical horizontal resolution of 5 km.Then, two intermediate models covering the shelf-slope domain and the coastal area surrounding the harbour are implemented with resolutions of 1 km and ~200 m, respectively.Finally, the harbour model is run with a resolution of 30 to 40 m, in order to solve the intricate harbour layout, which may include piers and docks (~30 m).The nesting ratio between different levels does not exceed 1:5 in order to avoid large resolution jumps that could filter out scales of interest (Debreu and Blayo 2008).
Each model provides initial and lateral boundary conditions for the model at the lower level with a prescribed frequency.In our case, the lateral boundary conditions are updated every hour.For harbours under tidal influence, this frequency is required to properly solve the tidal cycle.Conversely, in the Mediterranean some tests with different updating frequency have shown that updating at 6 and 12 hours leads to similar results inside the harbour (Grifoll 2009).Nevertheless, for simplicity we keep the 1-hour updating frequency for lateral boundary conditions.The preparation of the 3D fields for the model at the lower level (for initial and boundary conditions) is done as follows.First, we use a series of horizontal and vertical interpolations to transfer the outer solution variables needed for boundary and initial forcing to their respective positions on the inner grid.Second, a particular treatment is applied for points where coarse and high resolution bathymetries or coastlines do not match.In these cases, we compute the anomaly of the field with respect to the averaged value at each depth and then the anomaly value of the nearest point (in horizontal or vertical) is extrapolated to the problematic points.This procedure does not seem to introduce instabilities in the model and any errors in the procedure (i.e.interpolation/extrapolation errors) are expected to be significantly smaller than errors in the quality of the coarse resolution fields (see the System Validation section).This procedure is applied to the 2D and 3D fields.Additionally, a volume conservation constraint is applied to the current field.The velocity at the boundaries is corrected in order to ensure that the vertically averaged transport computed from the coarse resolution field and that computed from the interpolated fields have the same values (Penven et al. 2006).Finally, the atmospheric fields are also bilinearly interpolated on the ocean model grid and using only values over the sea.

Forcing fields
The atmospheric data for the three HCFS models were provided by the Spanish Meteorological Agency (AEMET, www.aemet.es).In particular, we used the 72-hour forecasts obtained from the High Resolution Local Area Model (HIRLAM, http://hirlam.org), with a spatial resolution of 0.16°.HIRLAM provides hourly fields of wind stress, sea surface temperature, surface net heat flux and surface freshwater (E-P).Wind observations were obtained from meteorological stations located in the harbour domain (see Fig. 1).River runoff was included in the regional applications, using data from the Global Runoff Data Centre (GRDC), which provides climatological runoff values for the most important rivers.In the near future, observed data for the Nervión river (Bilbao harbour) and Ebro river (Spanish Mediterranean area) will be used.For the Atlantic application (i.e.ESEOAT), 14 harmonic tidal components obtained from the MOG2D model (provided by LEGOS/POC) were introduced through the open boundaries.

Numerical codes
For historical reasons, there is no single model used at all the modelling levels.The system has had different development phases carried out by different teams, so different models have been used in the HCFS.However, it must be noted that all models are state-of-theart models that are comparable in complexity and performance.For clarity, we describe the numerical codes by regions (see Fig. 1 for model domains and Table 1 for a summary of model characteristics).

Mediterranean area
The forecasts of large-scale processes were provided by the Spanish Operational Ocean Forecasting System (ESEOO, Sotillo et al. 2007).In particular, for the Mediterranean area we used the ESEOMED product.ESEOMED is an implementation of the Die-CAST model (Dietrich/Centre for Air-Sea Technology, Dietrich et al. 1997), which covers the western Mediterranean Sea with 1/20° (~4-5 km) resolution.Two intermediate levels of increasing resolution were then implemented to reach the harbour scales.First, the northeastern Spanish shelf and slope area was covered with the SHE-CAT model at 1 km resolution.SHE-CAT was nested into ESEOMED with the aim of improving the representation of the mesoscale processes in the shelf-slope region.The dynamics in this area are dominated by a quasi-permanent slope current with strong variability (Millot 1999), and are affected by the presence of small eddies (Rubio et al. 2009), waves (Jordi et al. 2005) and filaments (Tintoré et al. 1990).Over the shelf, water circulation is mainly affected by wind-and density-induced currents, while tides have little influence (Bolaños et al. 2009).These processes have scales ranging from 1 to 20 km, so higher resolution than the one provided by the ESOMED model is required.Second, the shelf areas surrounding the harbours were modelled with a 200-m resolution coastal model (CST-BCN and CST-TGN for the Barcelona and Tarragona regions, respectively).The CST-BCN and CST-TGN models aim to improve the interactions of shelf processes with the coastline.Although SHE-CAT is able to reproduce the dominant processes, its resolution is too coarse to properly resolve the coastline details (capes or small bays).The model chosen for the Spanish Mediterranean shelf-slope area (SHE-CAT) and the coastal applications (CST-BCN and CST-TGN) was the SYMPHONIE model (Marsaleix et al. 2008).
Finally, the modelling of harbour dynamics requires a very high resolution implementation because of the usually complex harbour layout.In the Barcelona and Tarragona harbours, a 40-m mesh resolution was designed to correctly reproduce the harbour layout (PRT-BCN and PRT-TGN, respectively).The model chosen for these harbours was the Regional Ocean Modelling System (ROMS, Shchepetkin and McWilliams 2005).Some examples of ROMS implementations are given for a small-scale domain in Kim and Lim (2008), for coastal areas in Warner et al. (2010) and for harbour domains in Grifoll et al. (2009Grifoll et al. ( , 2011)).

Atlantic area
The ESEOAT application, another ESEOO operational product, was used to describe the large-scale processes in the Atlantic area.ESEOAT is an implementation of the Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS, Holt and James 2001).This application has a spatial resolution of 1/20° (~4-5 km) in the horizontal, and the vertical is discretized using sigma levels.The ESEOAT domain comprises the whole Atlantic Iberian shelf and slope, as well as the Bay of Biscay and part of the French continental shelf.The main characteristics of the water circulation within the Bay of Biscay are a weak oceanic circulation together with a frequent presence of eddies (Koutsikopoulous and Le Cann 1996).These slope water oceanic eddies (SWODDIES), with diameters of around 100 km (Pingree and Le Cann 1996), are the result of the continental margin current instabilities interacting with the bottom topography.Further, persistent poleward flow can be detected on the slope (Koutsikopoulous and Le Cann 1996).In the southern part of the Bay of Biscay, water circulation over the shelf is controlled mainly by wind-and density-induced currents (González et al. 2004), but because of the narrow width of the shelf (less than 40 km) a complex circulation pattern appears.Finally, density and tidal currents are also observed in the vicinity of the coastline, where frictional processes also gain importance (Grifoll et al. 2009).Obviously, ESEOAT is not able to reproduce the abovementioned processes near the coast because of its resolution, so two intermediate levels of increasing resolution were implemented: SHE-BIL (shelf-slope area) and CST-BIL (coastal area), with a spatial resolution of 1 km and 230 m, respectively.Finally, the Bilbao harbour circulation was resolved with a horizontal resolution of 43 m (PRT-BIL application).The ROMS model was chosen for SHE-BIL, CST-BIL and PRT-BIL.

INITIALIZATION STRATEGY
Initial conditions for each modelling level were provided by the model at the upper level (i.e.coarser resolution).This is the usual strategy in any nested system.It implicitly assumes that the coarse resolution model can provide a reasonable approximation to the actual field, so the initial conditions for the higher resolution model would be close to a stable state (i.e. in equilibrium with the model physics and forcings).This assumption is usually true in the open sea where higher resolution provides more details to the modelled fields but does not modify the general behaviour.However, for harbour modelling this is not the case.The complex geometry of the harbour is not reproduced by the coarser model.Harbour features have a typical size of a few tens of meters, while the coastal model has a resolution of 200 m.Therefore, the resolution provided by the coastal model to initialize the harbour model can be far from reality.It is therefore important to assess the impact on harbour forecasts of having wrong initial conditions.We must ascertain the typical time scale of the harbour system memory, which will determine the warming period (defined as the time needed by the model to "forget" the initial conditions) required by the harbour model before it can produce reliable forecasts.In this section, we present the sensitivity experiments run to address this issue.The analysis focuses on the Barcelona and Bilbao harbours as representative of a micro-and mesotidal environment, respectively.The strategy for analyzing the system memory was to run a set of simulations, each one starting from rest and forced by the same realistic conditions.The only difference between each simulation was the starting date.First, a reference experiment was run using realistic forcings: REF-BCN (October 2007) and REF-BIL (August 2008).Then, a set of three numerical experiments (perturbed simulations), with the initial date delayed, was carried out starting with null velocities and constant sea level.The delays used in the simulations were 5, 8 and 11 days for the Barcelona case and 1.5, 2.75 and 5.4 days for Bilbao.The differences between the delayed runs and the reference run reflect the impact of having wrong initial conditions on the model evolution.In order to quantify this, we computed the root mean square error (RMSE) between the delayed and the reference simulations every hour to obtain time series of model differences.The RMSE is defined as: where x 1 and x 2 are the velocity magnitudes at the n-points of the inner domain for the delayed and the reference simulations, respectively.The RMSE provides a useful indication of the forecast quality in an operational context.This diagnosis tells us if the delayed simulation converges to the reference simulation and, if it does, after how much time.In other words, it provides a measure of the system memory, given wrong initial conditions.The time evolution of the RMSE for the experiments carried out in the Barcelona and Bilbao harbours is shown in Figure 2. Additionally, the mean current is also plotted in order to have an idea of the relative importance of the errors.For the Barcelona case, starting from rest implies an initial RMSE of 2 to 3.5 cm s -1 , which is comparable to the averaged current in the harbour (relative errors ranging from 50% to 100%).There is an initial spin-up period of a few hours (between 6 and 12), during which the model first adjusts to an equilibrium solution.After this "numerical" spin-up, in which the difference with respect to the reference could be large, the RMSE decreases until it reaches a stable value of around 1 cm s -1 , after 72 to 96 hours.This value could be considered as a minimum forecast error induced by the initial conditions.The RMSE reduction is due to the effects of the forcings (wind and shelf currents), which dominate the model solution over the influence of the initial conditions.In Bilbao, the RMSE at the starting time is around 2-3 cm s -1 , about half of the averaged current.During the first 6 hours, there is a numerical spin-up to adjust to the model physics.Then, the RMSE decreases, reaching an almost constant value below 1 cm s -1 after only 48 hours.After 96 hours, the RMSE is reduced to 0.5 cm s -1 in all the cases.
It is interesting to note that the two harbours behave differently.Bilbao harbour has a shorter memory than Barcelona harbour, reaching the minimum RMSE faster.Also, this minimum RMSE is smaller for Bilbao harbour.The averaged currents and their variability are larger in Bilbao, so the relative impact of the errors is even smaller there than in Barcelona.Concerning the initialization strategy, the results suggest that the Barcelona harbour simulation should be launched 96 hours before the forecast period starts, in order to minimize the adverse effects of wrong initial conditions.In Bilbao, this warming period could be 48 hours.

SYSTEM VALIDATION
Once the system configuration has been fixed, it is important to compare the system results with observations in order to have an idea of the extent to which the forecasts (and therefore the derived products) are reliable.Unfortunately, sea observations were sparse in space and/or time, so an exhaustive validation of the models was not possible.However, we gathered all the information available during the system development phase and compared it with the model results whenever possible.Moreover, although the main focus of the HCFS is the water circulation in harbours, it is interesting to assess the impact of the quality of the regional and intermediate models on the quality of the harbour forecasts.Thus, we first validated the regional and intermediate models, and then the harbour models.

Atlantic area
The data available to validate the regional and intermediate models in the Atlantic area came from a drifting buoy deployed over the shelf in front of Bilbao.The buoy was kept at the sea surface between 3 and 6 June 2010.It started with a fast eastward displacement for 6 hours; it was then reversed for 8 hours and finally adopted a southeastward trajectory for 24 hours (see dotted line in Fig. 3a).In order to compare with this information, we analysed the trajectories of virtual buoys displaced by the surface current fields from the models.Both ESEOAT and SHE-BIL trajectories showed a similar southeastward displacement, although the total displacement of the SHE-BIL buoy was closer to the observations.The observed variability during the first 12 hours was not reproduced by any of the models, but the averaged behaviour was.It is interesting to note that, although the atmospheric forcing for ESEOAT and SHE-BIL was the same, the buoy trajectories were not, thus suggesting that the increase in model resolution helps to reproduce better the near submesoscale variability.
The Bilbao harbour results were validated using data from a measurement campaign carried out in summer 2008.In particular, we used velocity profiles from an Acoustic Doppler Currentmeter Profiler (ADCP) moored at the central river axis (see location in Fig. 1), from 23 July to 2 September 2008.The measured velocities were compared with the HFCS forecasts for the same period at the sea surface (0.5 m) and close to the bottom (7 m, see Fig. 4).The correlation coefficients for the surface layer velocities were 0.35 and 0.65 for the eastward and northward components, respectively.Also, it can be seen that the model underestimated the observed values (the slope of the fitted line is below 1).For the near-bottom circulation the model performance increased, showing correlation coefficients of 0.53 and 0.74 for the eastward and northward components, respectively.In this case, the model also underestimated the observed velocities.The averaged RMSE was 5.7 cm s -1 (see values in Fig. 4).It is worth mentioning that in the surface and bottom layers the largest velocities were in the northward direction (coherent with the tidal flow propagation), which is the component that was best solved by the model.The low correlations observed in the surface circulation are probably due to a wrong representation of the wind and the Nervión river runoff effects on the circulation, as pointed out by Grifoll et al. (2009), who also showed that the model properly solved the tidal induced circulation.Thus, errors in the wind fields or the river runoff could probably degrade the quality of the forecasts.The extent of this will be investigated in the next section.

Mediterranean area
The data used to validate the regional and intermediate models in the Mediterranean area came from the deep water observational network of Puertos del Estado (www.puertos.es)and corresponded to the area over the shelf in front of Tarragona (Fig. 1).Hourly observations of water velocity were obtained at 1 m depth from 1 to 31 October 2007.We computed the progressive vector from the observations and the models in order to summarize the velocity time series.The observations showed a clear southwestward direction due to the influence of the quasi-permanent northern current over the slope (Millot 1999).Superimposed on this general behaviour, there was some variability that could have been induced by slope current meandering, mesoscale structures or wind effects.The results from the ESEOMED regional model showed a weaker circulation oriented mainly westwards, thus differing from observations (Fig. 3b).This difference is due mainly to an erroneous representation of the slope current, which was too weak and/or even almost nil during part of the modelled period.Conversely, the results from the SHE-CAT intermediate model agreed better with the observations.The predominant velocity direction was closer to the observed one and the averaged velocity magnitude was the same as that observed.The increased spatial resolution in SHE-CAT, together with the improved vertical discretization, probably played an important role in the slope current representation.However, although the main behaviour was reasonably well-reproduced by the model, the variability was not.
The model was not able to reproduce at the correct time either the slope current meandering or the mesoscale structures (i.e.eddies and filaments).
The validation of the Barcelona harbour results involved considerable inconveniences.Measurements of water velocity inside the harbour were all obtained during an intensive field campaign in autumn 2003, while model results were not available until 2007.Thus, a direct comparison of observed and modelled results was not possible.However, we compared the statistical behaviour of the model in autumn 2009 with the statistical behaviour of observations in autumn 2003.Although particular events such as wind storms and mesoscale structures over the shelf were not the same in the two periods, we can expect these periods to show the same kind of variability on average.Therefore, we extracted the velocity time series from the model at the same locations where observations were obtained (see Fig. 1).The averaged magnitude of modelled currents inside the harbour was very similar to the observations (see Fig. 5).Near the north mouth (D1), modelled and observed currents were 7.1 and 7.5 cm s -1 , respectively.In the harbour channel (D2), the differences between model (12 cm s -1 ) and observations (8.9 cm s -1 ) were larger.Near the south mouth (D3), modelled (7 cm s -1 ) and observed (5.8 cm s -1 ) currents were closer.To compare the current variability in model and observations, we used the variance ellipses (Fig. 5).The size of the ellipses measures the magnitude of the variability; it can be seen that the model and observations showed similar ellipse sizes at all locations.Moreover, the model also correctly showed the polarization of current variability.At the north mouth, the current variability was oriented mainly in a perpendicular direction to the harbour mouth, as expected because of the harbour layout.However, it is interesting to note that variability in a parallel direction to the harbour mouth also had the right magnitude.In other words, the ellipse eccentricity matched both in the model and observations.In the harbour channel, the model ellipse was more elongated than in observations, as expected from the results shown above.This may be due to an incorrect representation of the harbour layout at this point.The channel width is narrow (only three model grid points), so small differences in this width imply significant differences in current magnitude.Finally, the south mouth variability was also correctly represented.Both modelled and observed ellipses were less polarized, reflecting a more isotropic variability.This is because the south mouth is considerably wider than the north mouth.In consequence, south mouth is less restrictive on current variability.In summary, we cannot say anything about the quality of the Barcelona harbour results for a particular period, but at least it seems that the model correctly reproduces the flow variability in terms of magnitude and variance.

SENSITIVITY OF HARBOUR FORECASTS TO FORCING ERRORS
In any operational system, the estimation of the forecast accuracy is almost as important as the forecast itself.Providing reliable estimates of the forecast errors would require an appropriate monitoring network, which is unfortunately not available at present in any harbours.However, what we can do is to perform sensitivity experiments to get a preliminary estimation of the forecast errors based on the known errors that are present in the forcing fields.In other words, we can use the difference between simulations as an estimation of the actual errors.In particular, here we focused on the potential impact of errors in the open boundary conditions and wind forcing.
The strategy was to compare pairs of simulations in which the difference between them is in the forcings (i.e.wind or lateral boundary conditions).In one case, we used the fields provided by the operational system (HIRLAM wind fields or lateral boundary conditions from the coastal model) and in the other measured values interpolated in the whole model grid.For the wind fields, we used the winds observed at the WS meteorological station (see Fig. 1 for its location) within the whole grid.For the lateral conditions we used the velocity profiles observed at all the boundary grid points and then we imposed a constraint to ensure volume conservation, as was done for the preparation of initial fields (see Nesting Strategy section).All the simulations were one month long.The Barcelona runs corresponded to October 2007 and the Bilbao runs to August 2008.Both periods involved a broad range of situations with alternate episodes of strong and weak winds and different shelf-current regimes (see Fig. 6a, b and Fig. 7a, b).Therefore, they can be considered representative enough of the variety of possible forc- ing situations, so the error statistics computed will be meaningful.Current and sea-level measurements outside Barcelona harbour were obtained from a coastal buoy of the XIOM network (www.xiom.cat) in October 2007, while those outside Bilbao harbour were obtained from an ADCP deployed in August 2008.The list of numerical experiments is summarized in Table 2 for the Barcelona (BCN) and Bilbao (BIL) harbours.In the reference runs BCN-1 and BIL-1, the harbour models were forced by observed winds and used the observed sea level and shelf currents for the lateral boundary conditions.In runs BCN-2 and BIL-2, the winds were provided by the HIRLAM forecasts while lateral boundary conditions were from observations.In runs BCN-3 and BIL-3, observed winds were used while lateral boundary conditions were provided by the coastal models.The difference between any pair of runs was quantified in terms of the RMSE at each  model grid point, computed for the whole simulation period.Additionally, we computed the noise-to-signal ratio (NSR), defined as: where σ is the standard deviation of the reference velocities.This quantity gives a measure of the relative importance of the errors compared to the field variability.Only surface velocities were used in the diagnostics for clarity and because this value is the most important field for harbour managers (i.e. it affects ship manoeuvring, oil spill dispersion, etc.).For these numerical experiments, the salinity and temperature in the water column were considered constants.

Sensitivity to errors in the wind field
Figure 6 (a, b) shows the observed and modelled (HIRLAM) wind time series.Wind variability in the two harbours is comparable in terms of wind strength (averaged wind stress of 0.05 Nm -2 ) and the number and strength of strong events (~5-6 events/month of 0.1-0.2Nm -2 ).The main difference is that strong wind events in Bilbao always correspond to eastward winds while in Barcelona they come from variable directions.Also, in Bilbao the sea breeze is not negligible, reaching up to 0.1 Nm -2 .The comparison between observed and modelled winds shows the same features in both harbours.Modelled winds reproduced most of the observed intense events although they systematically underestimated their strength.Errors in the modelled winds may reach up to 0.1 Nm -2 .Also, moderate variability was almost completely missed by the modelled winds.A likely explanation for the errors in the modelled winds would be that the atmospheric model spatial resolution (~14 km) was too coarse to consider topographic effects around the harbour domain (for instance, shelter provided by buildings).Also, winds were provided by mesoscale models, so perhaps they were not well suited for reproducing small-scale processes (~100 m) such as those affecting the harbour domain.
The differences between observed and modelled winds obviously translate to the modelled harbour circulation.The RMSE of water velocities (Fig. 6 (c,  d)) were similar in magnitude in both harbours, with an averaged value of ~1 cm s -1 and maximum values of 2.5 cm s -1 .The spatial structure of the errors was also comparable in both harbours.Errors were larger in the more exposed areas where the wind influence is greater.In Barcelona harbour, the influence of wind is in the central channel, where a continuous circulation between the two mouths is established.In Bilbao harbour, the wind influence is more noticeable near the mouth, which is also the widest part of the harbour.In both harbours, sheltered narrow areas are less affected by the wind.Even if they are under the same wind influence as the larger areas, the piers and docks inhibit the set-up of the wind-induced circulation.Also, in shallower areas (i.e.those close to land), friction is greater and wind-induced motions are damped.In consequence, in sheltered and shallow areas, the windinduced velocities were smaller, so the errors are small (below 0.5 cm s -1 ).
A complementary means for analysing the errors in the harbour circulation is the NSR (Fig. 6 (e, f)).In Barcelona harbour, the NSR due to errors in the wind field was between 40% and 50% in most of the domain.That is, errors were less than half of the variability.Near the south mouth, the NSR was below 20% while in some isolated areas it reached 60%.In Bilbao harbour, the NSR was highly anti-correlated with depth.Deeper areas showed a smaller NSR, with values below 20%, while shallower areas showed the highest NSR, reaching up to 100%.
Additionally, an analysis of the influence of the spatial wind resolution in the Barcelona harbour circulation is made.Wind fields provided by the operational meteorological models can be considered almost spatially homogeneous at the scales of interest.Their spatial resolution (5-10 km at best) cannot take into account the local topography around the harbour (i.e.buildings, hills, cliffs, etc.), so that the wind forecasts over the harbour have no spatial variations.A numerical experiment was run to check the impact of the lack of fine-scale structure on the wind fields (BCN-1b).In this case, we repeated the simulation forced by observed winds (BCN-1), but using spatially variable winds (see Table 2).These new wind fields were obtained by weighted averaging of wind data measured at meteorological stations (see Fig. 1) within the model grid.The weights were defined as a function of distance (w=exp(-D 2 /Lw), where D is the distance between the meteorological station and the model grid point and Lw=2 km).In both simulations (BCN-1 and BCN-1b), we used wind fields with the same averaged value but with a different spatial structure (homogeneous and spatially variable, respectively, see Fig. 7a,  b).The circulation in both experiments shows almost the same patterns inside the harbour and small differences between the two simulations are appreciated in the time evolution of the current magnitude (see Fig. 7c, d).The root mean square difference between the two simulations is 1.6 and 2.5 cm s -1 for the north and south mouths, respectively.The reason for such a good agreement is that harbour layout strongly conditions the preferential directions for the water motion, reducing the effects of the wind curl.

Sensitivity to errors in the lateral boundary conditions
The errors in the lateral boundary conditions in Barcelona harbour come mainly from the representation of the shelf currents.In Figure 8a, we show the currents measured in front of the harbour and the currents provided by the coastal model (CST-BCN).The magnitude, direction and temporal variability of observed currents are not reproduced.Observed currents show rapid changes in direction and maximum velocities of up to 20 cm s -1 .Conversely, the coastal model misses these variations in direction, and velocities are never over 10 cm s -1 .These large errors in the representation of shelf currents translate into large errors in the harbour forecasts (Fig. 8c), although they are located mainly near the south mouth and the central channel.This is consistent with the picture of the circulation induced by the shelf currents shown by Grifoll et al. (2011), who showed that southwestward shelf currents induce a circulation in the harbour from the north to the south mouth, with larger values in the narrowing of the channel.They also showed that the south mouth is strongly affected by current variability offshore, whatever its direction is.Maximum errors reach 7 cm s -1 in these areas, while in the rest of the harbour the values are lower than 1 cm s -1 (i.e.comparable to the errors induced by the wind).The NSR is again about 40% to 50% in a large part of the domain (Fig. 8e), but near the south mouth and in the central channel the NSR increases to 90%.Only in the most sheltered areas the NSR is below 10%.
In Bilbao harbour, the errors in the lateral boundary conditions are mostly related to the representation of the tides.Grifoll et al. (2009) showed that tides are the dominant driving mechanism of the Bilbao harbour hydrodynamics.Therefore, we compared the surface elevation measured outside the harbour with the values provided by the coastal model (CST-BIL).The coastal model provided acceptable results (see Fig. 8b), and discrepancies between the model and observations were in the amplitude (not the phase) and always less than 20 cm, that is, less than 10% of the total signal.However, these relatively small errors have a great impact on harbour circulation forecasts, with averaged RMSE ranging from 3 to 4 cm s -1 .The spatial structure of the RMSE is similar to that induced by errors in the wind fields (Fig. 8d).The largest values are located in the deeper areas and the smallest values are located in the shallower areas.This is consistent with the circulation pattern induced by tidal forcing (Grifoll et al. 2009), where larger velocities are found in the deeper areas of the harbour.Finally, the NSR distribution (Fig. 8f) is very similar to that shown in Figure 6, although the values are higher.Averaged NSR is 50% and values range from 30% in the deepest area to more than 100% in the shallower and sheltered areas.

DISCUSSION
The implementation of an operational system to provide harbour circulation forecasts is far from being a simple technical problem.The system requirements must be carefully assessed on the basis of the oceanographic characteristics of the environment and the particular conditions of the harbour design.For the Atlantic harbours, the tidal forcing is the main driving mechanism, while wind and shelf currents are of secondary importance (Grifoll et al. 2009).Conversely, for the Mediterranean harbours, where tidal influence is negligible, wind and shelf currents are the main driving mechanisms.However, it must be noted that in single-mouth harbours the shelf currents are much less influential (Mestres et al. 2007, Grifoll et al. 2011) and wind is the dominant driving forcing over short time scales (a few days).Additionally, in all harbours freshwater discharges may play a role, but they are considered to be a secondary forcing; in the case of a river inside the harbour domain (i.e.Bilbao harbour), its influence can be comparable to wind influence only during high runoff periods, although its main effect is on the thermohaline structure and it does not directly contribute to the circulation.In the case of a small-size river (i.e. that in Tarragona harbour) or the existence of other discharges such as wastewater, the impact is usually limited to a small area around the discharge point and only during the runoff period (Mestres et al. 2007, Grifoll et al. 2011).It is important to note that the relative importance of the driving mechanisms affecting harbour circulation may change from time to time.For instance, Grifoll et al. (2009) showed that under particular wind and river runoff conditions their contribution to surface circulation in Bilbao harbour is of the same order as tidal forcing.Also, Grifoll et al. (2011) showed that circulation in Barcelona harbour is strongly affected by southwestward shelf currents but almost unaltered by northeastward shelf currents.In consequence, all the elements that can potentially affect harbour circulation should be included in the system and their time variability should be correctly considered.
Wind forcing with reasonable spatial and temporal resolutions is routinely provided by meteorological services, so it is usually not a problem to include this forcing in the oceanic operational systems.Grifoll et al. (2011) showed that 3-hourly frequency for wind forcing is enough to reproduce the effects of wind on the Barcelona harbour circulation.Our experiments suggest that the spatial resolution of the wind field has little impact on the circulation inside the harbour.However, it must be mentioned that using higher resolution atmospheric models may lead to an improvement of surface winds in terms of both magnitude and temporal variability (i.e.reducing the errors observed in Fig. 6).The initial and lateral boundary conditions are more difficult to obtain because there are few opensea operational systems, and when they exist their spatial resolution is far from optimal for properly solving coastal circulation.Therefore, if currents or sea level over the shelf must be included in the system (i.e. at Bilbao and Barcelona harbours), a downscaling of open sea operational products must be performed.Our choice was to use different nesting levels to transfer the information from the open sea to the coastal domain, including the particular features of the coastal area.For Tarragona harbour, this would probably not be necessary, provided that shelf currents have little influence on the inner harbour hydrodynamics (Mestres et al. 2007).However, because the intermediate models are already available for this area, we also included the shelf current forcing in Tarragona.Finally, river discharges are difficult to forecast and at present only climatological values are used.
Once the system requirements have been assessed and the system is implemented with all the required elements, it must be configured.First, the model parameters must be calibrated in order to compare model results with observations.In our case, because limited data were available, this was not possible and typical values used in other regions were used.Second, the operational functioning must be established.In other words, we have to decide the forecast horizon and the restart procedure.The forecast horizon is determined by the availability of the forcing fields, which is 72 hours for the HIRLAM wind fields and the ESEOO regional models (ESEOMED and ESEOAT).The restart procedure is closely linked to the system memory.We have seen that the impact of having wrong initial conditions is noticeable until 96 hours after the start in the Barcelona case, and 48 hours in the Bilbao case.In Tarragona harbour, we made no estimation, but as it is a small harbour driven only by the wind, it is probably a safe choice to use the same results as those obtained for Barcelona.Therefore, the harbour models are configured to start 96 hours (48 hours in the case of Atlantic harbours) before the forecast phase begins.
Although this study has been carried out on three harbours, we have enough information to be able to make some recommendations for implementing a similar system in other harbours.Flow patterns inside semienclosed domains with intrincate shorelines are usually complex (e.g.Cucco and Umgiesser 2006, Niedda and Greppi 2007, Ribbe et al. 2008).Moreover, studies on harbour circulation highlight the 3D behaviour of the hydrodynamics inside the harbour and the interaction with winds, shelf currents and freshwater sources.Therefore, the first recommendation is that a 3D primitive equation model with enough spatial resolution to correctly represent the complex harbour layouts must be used.The second recommendation concerns the forcing fields.In Atlantic harbours (i.e.harbours in a mesotidal environment), wind forecasts and a good representation of tides outside the harbour (i.e.sea level and currents) must be provided to the system.Tidal forcing could be provided by a coastal 3D model (as in our case) or simply by a tidal model, which is simpler and usually accurate enough.In this case, effects of non-tidal shelf currents would not be included in the system, but such effects are considered to be of secondary importance.In Mediterranean harbours (i.e.harbours in a microtidal environment), wind forecasts are required and if the harbour has two mouths, shelf currents should also be provided.In all cases, having river runoff forecasts (if there is a river inside the harbour domain) would improve the system quality.However, because they are difficult to obtain, climatological values or past observations could also be used, assuming that rivers are less influential on harbour circulation.Small rivers (i.e. the Francolí River in Tarragona harbour, with an averaged runoff of 2 m 3 s -1 ) have a small impact on circulation (Mestres et al. 2007), whereas larger rivers (i.e. the Nervión River in Bilbao harbour, with an averaged runoff of 35 m 3 s -1 ) can moderately affect the circulation but only close to the river mouth (Grifoll et al. 2009).
A proper assessment of the quality of system results would require an extensive dataset which is not available at present.However, with the limited data available we can draw some conclusions about the HCFS quality.The intermediate models (SHE-CAT and SHE-BIL) provide an improvement of the regional models (ESEOMED and ESEOAT).The mean circulation improves in both magnitude and direction (see the validation section), but it is still not able to reproduce the high-frequency variability (order of days).For the harbour models, we were able to compare the model results with contemporary data (i.e. from the coastal model to force the harbour model and observations within the harbour) only in Bilbao.The averaged RMSE is about 6 cm s -1 and the averaged correlation is ~0.6, which is reasonably good considering the system complexity.In Barcelona harbour, we were only able to perform a statistical analysis and it seems that the model correctly reproduces the typical circulation patterns inside the harbour.Unfortunately we were unable to go further in the Barcelona model validation.In Tarragona harbour, Mestres et al. (2007) showed that the general behaviour of the model matches current observations, but errors may be large during certain periods and errors in the current field may vary between 2 and 15 cm s -1 , depending on wind conditions.In summary, the HCFS can provide a reasonable approximation to reality, particularly to the averaged behaviour, but the range of uncertainties is still large.
The forecast errors may come from different sources than the initial and boundary conditions or atmospheric forcing.The sensitivity experiments carried out here give us an idea of the relative importance of different error sources.Having wrong initial conditions induces large errors in the initial simulation period.However, after the warming period (48 hours for Atlantic harbours and 96 hours for Mediterranean harbours), the RMSE averaged in the whole domain is stabilized at around 1 cm s -1 .The errors in the wind forcing have more influence in areas where wind-induced circulation is stronger, that is, in the wider areas such as near the harbour mouths or in large basins.In small docks, the circulation is much weaker because the harbour layout prevents the establishment of the wind-induced circulation, so the magnitude of errors is less important there.From our experiments, typical RMSE values due to wind inaccuracies are around 3 to 4 cm s -1 .Finally, with expected RMSE values of 4 to 7 cm s -1 , the errors in the information prescribed through the lateral boundary conditions are the most significant in Barcelona and Bilbao harbours (in Tarragona, lateral boundary conditions do not affect harbour circulation).This finding has several causes.First, lateral forcing is the dominant driving mechanism in Bilbao (Grifoll et al. 2009) and Barcelona harbours (Grifoll et al. 2011), so the system is more sensitive to errors in these fields.Moreover, in Barcelona harbour, errors in the lateral forcing (i.e. shelf currents, see Fig. 8a) are especially large.Additionally, model parameters or limited phys-ics in the model may have an influence on harbour forecast quality.This influence has not been analyzed in this paper but some sensitivity studies carried out by Grifoll (2009), varying different model parameters (such as bottom friction, and horizontal and vertical viscosity), suggest that their impact on harbour circulation is secondary when compared to the impact of errors in the wind field or lateral boundary conditions.An improved error characterization may be achieved when monitoring networks are deployed in the harbours.Alternatively, advanced techniques such as stochastic modelling (see, for instance, Jordà and De Mey [2010] for an application in the western Mediterranean) could be used to improve error characterization.
In summary, we present here the first operational forecasting system for harbour circulation developed for Spanish harbours.To our knowledge, no other similar initiatives have been developed in southern Europe.At present, the system is far from perfect, but for the first time a tool related to water circulation can be used by harbour managers and may be of crucial importance for ship manoeuvring operations and oil spill control.Furthermore, circulation forecasts are the basis of several methodologies for water quality degradation management (Marin et al. 2008, Choi et al. 2009, Grifoll et al. 2010).Guerra-García et al. (2005) also showed that sediment quality and benthonic biodiversity in a harbour are closely related to harbour circulation.Some authors have suggested that in nested systems the quality of shelf models strongly depends on the quality of the information provided in its open boundaries (Auclair et al. 2001, Barth et al. 2008, Mason et al. 2010).In our case, simulation length is probably too short to notice the same problems that those authors have reported.Nevertheless, our sensitivity experiments show that an improvement in the quality of the shelf models will directly improve the quality of the harbour model.Therefore, efforts should be devoted to improving the forecast quality at the different nesting levels and not only at the harbour level.Future developments of the system should focus on the improvement of forcing fields (e.g.improving the quality of atmospheric fluxes or the operational systems that provide lateral boundary conditions), the implementation of monitoring networks (e.g. for validation and calibration tasks), and the improvement of the modelling system (e.g. through the use of two-way nesting or data assimilation).

CONCLUSIONS
The Harbour Circulation Forecasting System developed for Spanish harbours is presented in this paper.A preliminary analysis of the driving forcing mechanisms in different harbours has been made in order to determine the system requirements.A four-level nested system of ocean models, extending from the regional scale to the harbour scale, has been implemented.Sensitivity experiments oriented towards determining the optimal initialization strategy suggest that, for the Atlantic harbours, model runs should start at least 48 hours before the forecast phase.For the Mediterranean harbours, the warming phase should last 96 hours.The system validation shows encouraging results for the harbour circulation forecasts.In general, the models reproduce the main circulation patterns, although they miss particular features.This is due mainly to the limitations of the forcings.Additional sensitivity tests have shown that errors in the lateral boundary conditions are the most influential in the quality of the harbour forecasts.The reason is twofold.First, the conditions over the shelf (currents and elevation) are the dominant forcing in harbours in tidal environments or with two mouths.Second, the quality of the information on those conditions is significantly lower than the quality of the wind fields.Nevertheless, errors in the wind fields induce errors in the harbour forecasts of about 30%, so they are not negligible either.Although the system is far from perfect, this is the first time that harbour managers have had access to operational products based on circulation predictions, which is of crucial importance in terms of ship operations, risk management and oil spill control.

Fig. 1 .
Fig. 1. -Top: Barcelona, Tarragona and Bilbao harbour locations (PRT-BCN, PRT-TGN and PRT-BIL, respectively) on the Spanish coast and the domains of the intermediate models (i.e.SHE-CAT, SHE-BIL, CST-BCN, CST-TGN and CST-BIL).The star indicates the location of the measurement point used to validate the SHE-CAT shelf model.Bottom: harbour layouts.The circles and crosses in Barcelona and Bilbao harbours indicate the location of the ADCP and wind measurements, respectively.WS is the reference wind station for Barcelona harbour.

Fig. 2 .
Fig. 2. -Time evolution for the RMSE (cm s -1 ) between the delayed and the reference (REF-BCN and REF-BIL) simulations for runs started on different dates.(a) Results for Barcelona harbour; dotted, black and grey lines correspond to runs started on 5, 8 and 11 October 2007, respectively.(b) Results for Bilbao harbour; dotted, black and grey lines correspond to runs started approx.on 1, 3 and 5 August 2008, respectively.The time evolution of the mean current(cm s -1 ) within the harbours (black thin lines) is also shown.

Fig. 3 .
Fig. 3. -(a) Observed tracking of a surface drifting buoy near Bilbao harbour (dotted line) compared with the predicted trajectories by ES-EOAT (black line) and SHE-BIL (grey line) models, between 3 and 6 June 2010.(b) Progressive vectors from the ADCP (grey line) located over the Mediterranean shelf and from ESEOMED (black line) and SHE-CAT (dotted line) models, between 1 and 31 October 2007.

Fig. 4 .
Fig. 4. -Comparison between the measured and computed U and V current components at the sea surface (a, c) and 7 m depth (b, d), at the location of the ADCP in the Bilbao harbour, for the period between 23 July and 2 September 2008.The correlation coefficients (R) and the RMSE (cm s -1 ) for each component are also shown.

Fig. 5 .
Fig. 5. -Variance ellipses of measured (in black) and modelled (in grey) currents for the autumn period of 2003 and 2009, respectively, in Barcelona harbour.Arrows show mean velocities.D1, D2 and D3 denote the location of the measurement points referred to in the text.

Fig. 6 .
Fig. 6. -Sensitivity of circulation to wind errors in Barcelona (left panel) and Bilbao (right panel) harbours, for October 2007 and August 2008, respectively.(a, b) Wind stress components (Nm -2 ) used for the experiments: observed winds (blue line) and modelled winds obtained with HIRLAM (red line).(c, d) RMSE (cm s -1 ) obtained with HIRLAM winds.(e, f) NSR (%; see text for details).The results are masked outside the harbour because the colour scales are too different (larger outside).

Fig. 7 .
Fig. 7. -Snapshot of the wind forcing used in the experiments BCN-1 (a) and BCN-1b (b), corresponding to 3 October 2007.Time series of current magnitude in the north (c) and south (d) mouths for both experiments.A colour version of this figure may be found in the online electronic manuscript.

Fig. 8 .
Fig. 8. -Impact of lateral boundary conditions errors in Barcelona (left panel) and Bilbao (right panel) harbours, for October 2007 and August 2008, respectively.(a) Shelf current components (cm s -1 ) from observations (blue line) and provided by the coastal model (red line).(b) Sea level (m) from observations (blue line) and provided by the coastal model (red line).Also, the difference is plotted (green line).(c, d) RMSE (cm s -1 ) obtained with the coastal model results for the lateral boundary conditions.(e, f) NSR (%; see text for details).The results are masked outside the harbour because the colour scales are too different (larger outside).

table 1 .
-Characteristics of the numerical models implemented.

table 2 .
-Numerical experiments designed to assess the impact of the forcing errors on the harbour forecasts.