ROMS Based Hydrodynamic Modelling Focusing on the Belgian Part of the Southern North Sea

Abstract

The North Sea is a shallow sea that forms a complex physical system. The nonlinear interaction of the astronomical tides, varying wind fields and varying pressure systems requires appropriate approaches to be described accurately. An application based on the advanced numerical model Regional Ocean Modeling System (ROMS) was newly developed by the authors, tailored to simulate these hydrodynamic processes in the North Sea and the Belgian Continental Shelf, which is the area of particular interest in the present study. The purpose of this work is to develop and validate a state-of-the-art three-dimensional numerical model to form the basis of a compound operational and forecasting tool for the Belgian coastal zone. The model was validated with respect to water levels and temperature. Validation for astronomical tides was accomplished through the comparison of the principal constituents between the model results and observations at a number of tidal gauges in Belgium and other countries. A statistical analysis of the results showed that the model behaves as expected throughout the North Sea. The model response to the varying meteorological conditions was also validated using hindcast data for 2011 as input. In this case, the comparison between observed and modelled water levels showed a good agreement with average RMSE in Belgium 9.5 cm. Overall, the added value of this work is the development of an independent model for validation and comparison with other models and which can be used as an efficient tool for operational and forecasting purposes.

1. Introduction

The North Sea spans an area of more than 970 km in latitude and 580 km in longitude and is surrounded by eight highly developed industrial countries. The life of more than 50 million people is directly affected by the presence of this marginal sea, while more than 200 million live in the catchment area of the rivers discharging into it. Due to its importance to the highly industrialized area, the North Sea is one of the most intensely investigated seas in the world [1,2]. However, the strongly nonlinear character of the processes governing the North Sea and the existence of measured data gaps referring to the spatial distribution of significant parameters (water velocities, salinity, temperature, precipitation) stimulate continuously the interest of the scientific and engineering community for a more in-depth understanding of its complex environment.

Large databases with reference to the three-dimensional (3D) and temporal distribution of oceanographic and meteorological observations are maintained over the decades by major data centers for the North Sea, e.g., the British Oceanographic Data Centre (BODC), the German Oceanographic Data Centre (DOD), the Service Hydrographique et Océanographique de la Marine (SHOM), the Data Publisher for Earth & Environmental Science (PANGAEA), etc. These databases are constantly enriched with new in situ observations and remote sensing information, mainly provided by governmental agencies in the states surrounding the North Sea, e.g., the Ministry of Infrastructure and Water Management (Rijkswaterstaat) in the Netherlands, Federal Maritime and Hydrographic Agency of Germany (BSH), Danish Meteorological Institute (DMI), etc.

Nevertheless, the observational datasets alone are not sufficient for forecasting applications. Thus, a number of hydrodynamic numerical models have been developed and applied over the North Sea domain for both research and operational purposes. Namely well-established models are DELFT3D [3], MIKE 3 [4], DCSM [5], TELEMAC-3D [6], HAMSOM [7], NORWECOM [8], POLCOMS [9], COHERENS [10], HBM [11], TRIWAQ [12], MARS 3D [13] and GHER 3D [14]. Some of these models are applied on an operational basis for forecast services and used for predicting storm surges and tidal elevations (BSH, NORWECOM). Other models are used as consulting tools in the process of decision making, e.g., MIKE3, DELFT3D, etc.

In particular, the aforementioned models have been applied and proven efficient to reproduce the main features of the North Sea circulation [15]. However, differences exist between the various hydrodynamic models that result from the various physical, mathematical and numerical formulations that they embed. The evaluation of the different models behavior in the North Sea has been accomplished through intercomparisons performed in the context of various projects and working groups [15,16,17,18,19,20]. The general conclusion from these studies was that it is extremely difficult to consistently quantify and separate the causes of the differences in the model results. No single model parameter (spatial resolution, turbulence closure scheme, model domain, etc.) can explain the different model results, nor is the RMSE a reliable measure of such comparisons if long time periods are not considered, i.e., at least one calendar year [5,21]. The various hydrodynamic models are developed for different purposes, leading to different formulations and different implementations [22].

As ships become bigger and, therefore, ship drafts become deeper, in addition to the limited amount of waterways in the Belgian part of the North Sea, a new challenge is formed. This refers to the optimization of the use of the available waterways, while maintaining safe ship navigation and minimizing the amount of expensive dredging campaigns. Next to safety of navigation, optimizing coastal protection also benefits from the seamless forecasting of the water levels. Furthermore, the enhanced climate change phenomena intensify the variability of sea level and coastal geomorphology. In the present study, we tailored the Regional Ocean Modeling System (ROMS) to simulate hydrodynamics due to astronomical tides and meteorological conditions in the North Sea and, particularly, in the Belgian Continental Shelf (BCS). ROMS is a state-of-the-art 3D model and, to the best of the authors’ knowledge, there are very few research works presenting its application in the North Sea. These works have the particular interest on estimating the tidal-stream renewable energy and they used validation locations only inside the United Kingdom [23,24]. There are also studies referring to modelling eco-hydrodynamics within the BCS with other models, e.g., COHERENS, TELEMAC-2D, NEVLA, GHER [25,26,27,28].

The motivation of the present research was to develop a new hydrodynamic model with particular interest in optimizing its accuracy in the Belgian Continental Shelf, for research or operational purposes. The added value is threefold. At first, the independent model development offers a tool for validation of other models applied in the North Sea and the Belgian Coastal Shelf. Secondly, the wide range of capabilities, options and modules included in ROMS offer the opportunity to use the model for more focused research on the accuracy of the input, processes and output. Thirdly, providing a validated hydrodynamic ROMS model has the merit to form a basis for further study using a research type model. The present model is validated with respect to tidal constituents, water temperature and statistical errors using hindcast atmospheric data for 2011. Its behavior is evaluated in the light of forming the basis of a compound tool for forecasting and operational purposes, mainly within the Belgian coastal zone. ROMS has the advantage of being an integrated system that also includes modules for sediment transport and morphological evolution. Hence, the developed numerical tool will be further extended in the future to account for coastal sedimentation and beach erosion in the Belgian coastal zone.

2. Dynamics of the North Sea

The North Sea is a shallow sea adjacent to the North Atlantic with a mean depth of approximately 80 m. The bathymetry varies from 200 m at the north entrance to approximately 50 m over an area extending from the Dogger Banks to north Denmark. The depth values are further reduced to less than 20 m along the Dutch-German coast, while the highest values of about 800 m are met in the Norwegian Trench. To the west, the North Sea connects to the Atlantic Ocean through the English Channel, which is maximum 100–120 m deep along its centerline. At the Strait of Dover, between Dover and Calais, the bottom is about 45 m deep. A map of the domain of interest is shown in Figure 1.

The broad connection of the North Sea with the surrounding ocean, as well as the intense continental impacts of NW Europe, establish a complex physical regime that requires appropriate modelling approaches to be described accurately. In particular, the governing processes result from the nonlinear interaction of various factors, i.e., the astronomical tides, the North Atlantic Oscillation (NAO), the low pressure systems, the varying wind fields, the freshwater discharges, the continental heat flow, with the complex coastal geography and the relatively shallow waters. The NAO is a phenomenon in the North Atlantic that includes oscillations of the high and low pressure systems in low and high geographical latitudes, respectively. These decadal oscillations are transferred to the North Sea, mainly, through the atmosphere and the control of the westerly winds, and less through direct exchange of water masses [29].

The tidal motion in the North Sea is not the result of direct tide-generating forces, but the co-oscillating response of tidal waves generated in the North Atlantic Ocean [30]. This motion is principally semi-diurnal and progresses cyclonically around the North Sea, with the largest amplitudes along the English Channel, the coasts of eastern England and the German Bight. The main tidal constituents enter the North Sea through an open boundary with the North Atlantic Ocean in the north and the Dover Strait in the south. The latter has a length of approximately 35 km and may function as a source or sink of tidal energy, which may cause a phase shift of harmonic constituents. On the other hand, the Skagerrak Strait, which connects Norway and Sweden with the Jutland peninsula of Denmark, is approximately 100 km wide and mainly functions as a sink of tidal energy [1]. The most dominant tidal constituent in the North Sea is M2-tide, which propagates in a counterclockwise direction through the main basin. Near the coast of Scotland and England it propagates almost as a perfect Kelvin wave. The spatial pattern of the M2-tide is characterized by three elevation amphidromic points in the North Sea. One is located in the Southern Bight, a second one northeast of it, in the Dogger Banks, and the third one even further to the northeast, near the Norwegian coast.

Through the vertical flux of momentum the atmosphere significantly controls the general circulation of the North Sea. Variations in atmospheric pressure and especially strong winds force currents and associated elevations in the extensive shallow waters of the North Sea. In central-northeastern areas, these may be larger than tidal motion, and comparable elsewhere on occasion. The prevailing winds on the northwest European Shelf are westerly which cause an intense cyclonic circulation, while for easterly winds this pattern reverses. Low pressure systems are also related to storm surges in the North Sea and extreme events [29]. In order to simulate the hydrodynamic behavior of storm surges the input from high spatial resolution meteorological models is required.

The Belgian coastal zone is approximately 65 km long and extends 87 km offshore. The average and maximum depths are 20 m and 45 m, respectively. In the cross-shore direction the seabed is mildly sloping, ranging between 1:20 and 1:90 [31]. However, the bathymetry is also marked by a complex system of large and elongated submerged sandbanks, the Flemish Banks. Their orientation is, in general, alongshore and their length may be up to 15 km. The crest to trough depth differences may be up to 20 m [32]. The meteorological conditions over the Belgian Continental Shelf (BCS) are highly correlated with the North Atlantic Oscillation (NAO) [33] performed a simple assessment of the relative importance of the forcing from wind stresses and atmospheric pressure gradients. In their analysis they neglected the variations of water depth and the possibility of surges generated by atmospheric pressure gradients in deep water, which may propagate into the shallow sea. By additionally considering typical values for the involved parameters, they concluded that storm surges in this region are wind-dominated. An annual cycle with generally lower winds during summer is identified over the Belgian coastal zone. In addition to this, a significant variability of both wind speed and direction is observed at time scales of hours to days. This is mainly associated with the passage of low pressure atmospheric systems.

The dynamics of the Belgian coastal zone is determined at a large extent by the strong semi-diurnal tidal signal. For example at Ostend the amplitude of the M2 and S2 tides reach 1.80 m and 0.60 m, respectively. Furthermore, the shallow topography of the Belgian Continental Shelf enhances the nonlinear effects due to bottom friction and tidal-induced changing depths. These processes are more intense than nonlinearity due to advection. Consequently, tidal curves are severely non-harmonic as significant energy is transferred to higher-order harmonics, e.g., M4, M6, MS4, etc. [29].

As a result of the large tidal currents, the water column remains relatively well-mixed all through the year. In particular, no permanent haline stratification can be associated with the plume of the Scheldt River, because of the low freshwater discharge and the strong tidal mixing. Similarly, no thermocline is formed in the BCS due to the intensity of tidal mixing.

3. Model Configuration

3.1. Mathematical Background

The Regional Ocean Modeling System (ROMS) is a three-dimensional, free-surface, terrain-following numerical model that solves the Reynolds-Averaged Navier-Stokes equations using the hydrostatic and Boussinesq assumptions [34,35,36,37,38,39]. The mathematical formulation of ROMS in Cartesian form consists of the continuity equation,

$$\nabla ·\mathit{v}=0,$$

(1)

the horizontal momentum equations,

$$\frac{\partial u}{\partial x}+\left(\mathit{v}·\nabla \right)u-fv=-\frac{\partial \phi }{\partial x}-\frac{\partial }{\partial z}\left(\overline{u^{\prime }{w}^{\prime }}-\nu \frac{\partial u}{\partial z}\right)+F_u+D_u,$$

(2)

$$\frac{\partial v}{\partial y}+\left(\mathit{v}·\nabla \right)v+fu=-\frac{\partial \phi }{\partial y}-\frac{\partial }{\partial z}\left(\overline{v^{\prime }{w}^{\prime }}-\nu \frac{\partial v}{\partial z}\right)+F_v+D_v,$$

(3)

the vertical momentum equation,

$$\frac{\partial \phi }{\partial z}=-\frac{\rho \mathrm{g}}{{\mathsf{\rho }}_{\mathrm{o}}},$$

(4)

the scalar transport equation for salinity, $S$, and potential temperature, $T$,

$$\frac{\partial c}{\partial t}+\mathit{v}·\nabla T=-\frac{\partial }{\partial z}\left(\overline{c^{\prime }{w}^{\prime }}-{\nu }_{\theta }\frac{\partial c}{\partial z}\right)+F_c+D_c,$$

(5)

and the equation of state,

$$\rho =\rho \left(T,S,P\right)$$

(6)

where the 3D flow velocity has components $\mathit{v}=\left(u,v,w\right)$ in the $x$, $y$, and $z$ directions, respectively, $\rho$ is the water density, $P$ is the total pressure, $t$ is time, $c=\left\{T,S\right\}$, $\nabla =\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y}\right)$ is the horizontal gradient operator, $f$ is the Coriolis parameter, $\phi =P/{\mathsf{\rho }}_{\mathrm{o}}$ is the dynamic pressure, ${\mathsf{\rho }}_{\mathrm{o}}$ is a reference density, $\nu$ is the molecular viscosity, $F_u,F_v,F_c$ are forcing/source terms, $D_u,D_v,D_c$ are horizontal diffusive terms resulting from Reynolds averaging, $\mathrm{g}$ is the gravitational acceleration, ${\nu }_{\theta }$ is the molecular diffusivity, primes (′) refer to turbulent variations of the relevant variables, and overbar (^―) denotes time-averaged values.

The above equations are closed by parametrizing the Reynolds stresses and turbulent tracer fluxes as:

$$\begin{gathered} \overline{u^{\prime }{w}^{\prime }}=-K_M\frac{\partial u}{\partial z},\overline{v^{\prime }{w}^{\prime }}=-K_M\frac{\partial v}{\partial z},\overline{c^{\prime }{w}^{\prime }}=-K_c\frac{\partial c}{\partial z}, \\ \overline{u^{\prime }{v}^{\prime }}=-A_M\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right),\overline{c^{\prime }{u}^{\prime }}=-K_{c,h}\frac{\partial c}{\partial x},\overline{c^{\prime }{v}^{\prime }}=-K_{c,h}\frac{\partial c}{\partial y} \end{gathered}$$

(7)

where $K_M$ and $A_M$ are the vertical and horizontal eddy viscosity for momentum, respectively, and $K_c$ and $K_{c,h}$ are the vertical and horizontal eddy diffusivity for tracers, respectively.

Equations (1)–(6) are expressed in Cartesian horizontal coordinates. However, in the present application spherical coordinates were used, which implied the presence of additional metric terms [35]. The model state variables are staggered using an Arakawa C-grid. In the vertical plane, a generalized vertical, terrain-following, $S-$ coordinate system is used which essentially flattens the variable bottom at $z=-h\left(x,y\right)$, where $h\left(x,y\right)$ is the still water depth. The nonlinear vertical transformation functional, $S$, is given by:

$$S\left(x,y,\sigma \right)=\frac{h_{\mathrm{c}}\sigma +h\left(x,y\right)\mathrm{C}\left(\sigma \right)}{h_{\mathrm{c}}+h\left(x,y\right)},$$

(8)

where $\sigma$ is the vertical distance from the surface measured as the fraction of the local water depth, $\mathrm{C}\left(\sigma \right)$ is a nondimensional, monotonic, vertical stretching function ranging from $-1\le \mathrm{C}\left(\sigma \right)\le 0$, and $h_{\mathrm{c}}$ is a positive thickness controlling the stretching [40,41,42,43].

In the horizontal plane, a centered, second-order finite-difference approximation is adopted. Time integration is performed using a split-explicit time step, integrating the depth-integrated equations with a shorter time-step (barotropic mode) than the full 3D equations (baroclinic mode). Time-stepping is accomplished through a LeapFrog-3rd order Adams-Moulton (AM3) predictor-corrector scheme. The algorithms that comprise the ROMS computational kernel are described in detail in [41,44,45].

3.2. Study Area—Computational Grid

The study area covers the entire North Sea bounded by the English Channel to the west, the Atlantic Ocean and the Norwegian Sea to the north and the Skagerrak Strait to the northeast. The Belgian Continental Shelf is of particular interest in this study and, thus, a refined grid (BSC Model) is nested within the larger North Sea model (NS model). The latter spans the area between 48.5° N and 59.5° N and from 4.5° W to 9° E. The NS model domain was discretized using a horizontal curvilinear grid, with longitudinal and latitudinal resolution of 1/12° and 1/22°, respectively, resulting to a grid resolution of approximately 5.0 km × 5.0 km. This resolution was selected as a balance between accuracy, numerical stability and computational efficiency. An ideal 3D model of the North Sea should have horizontal and vertical grid sizes of less than about 5 km and 10 m, respectively, to be able to resolve the horizontal and vertical density gradients and the depth varying residual current patterns. Moreover, in order to resolve the Rossby radius of deformation throughout the considered domain, a grid size of maximum 5 km is required. A grid size of 2.5 km × 2.5 km was also investigated but the model results were similar in the greatest part of the domain (outside the Belgian and Dutch Continental Shelf), while the computational time was almost four times higher. On the other hand, the below described refinement was necessary for accurately predicting the nonlinear tidal interactions and atmospheric effects within the shallow BCS. The model bathymetry was based on the ETOPO [46] global bathymetric dataset, which is available at a resolution of 1 arc-minute (see Figure 2).

The BCS Model covers the Belgian Continental Shelf, the southeast part of United Kingdom, the northeast part of French coast, and the south of Netherlands, just north of Scheldt River estuary, bounded between 50.8° N and 52.0° N and from 0.7° E to 4.4° E. The longitudinal resolution is 1/63°, while the latitudinal is 1/110°. Thus, the refinement ratio between parent and child grid is 1:5, resulting to a horizontal grid size of approximately 1.0 km × 1.0 km for the nested model. A fully online nesting procedure is applied for the two models. The bathymetric data for the BCS is a collection of bathymetric survey datasets and were merged with the ETOPO database. Terrain Ruggedness Index was employed to check the resulting merged bathymetry on internal consistency [47] (see Figure 3). The outline of the domains of the two models is shown in Figure 1.

For both NS and BCS models the fully 3D (baroclinic) formulation was developed since the resulting computational load was not prohibitive, i.e., 2.5 h and 41 h, respectively, for a full year simulation on a 72-threaded server. In particular, 15 sigma vertical layers were applied throughout the computational domain. The vertical stretching parameters involved in $\mathrm{C}\left(\sigma \right)$ were set ${\theta }_s=8.0$ and ${\theta }_b=4.0$ [42].

3.3. Model Settings

Southern North Sea is, in general, a well-mixed environment and only locally, e.g., Scheldt estuary, increased mixing of fresh and sea water is taking place. Hence, in the present study the sea water temperature, salinity and density were assumed uniform and time invariant, $T={12}^{\circ }\mathrm{C}$, $S=36‰$ and $\rho =1025$ kg/m³, respectively. For the bottom friction the conventional quadratic formulation was employed with a constant friction coefficient $\gamma =0.004$. The Coriolis effect was also taken into account. As far as the air-sea interaction is concerned, the surface boundary layer was described based on the bulk parameterization by [48]. The effect of clouds and rain were ignored, while the surface air humidity and temperature were assumed invariant and uniform. Furthermore, among all rivers discharging into the North Sea, only Scheldt was taken into account since it is the only river near the area of particular interest, i.e., the Belgian Continental Shelf.

The applied vertical mixing parameterization relied on the Generic Length Scale (GLS) approach for turbulence closure [49]. Among the various formulations of this two-equations model, the one that corresponds to the traditional $k-\epsilon$ model was implemented, in combination with the stability function described by [38,50]. In the horizontal plane, a viscosity coefficient with a maximum value of 5.0 m²/s was considered for momentum diffusion.

At the open boundaries the implicit Chapman condition was imposed for the free surface elevation [51]:

$$\frac{\partial \zeta }{\partial t}=\pm \sqrt{\mathrm{g}h}\frac{\partial \zeta }{\partial \xi }$$

(9)

where $\zeta$ is the free surface elevation and $\xi$ is the curvilinear coordinate.

With regard to the barotropic velocities a Flather boundary condition was applied [52]. For the normal component of the barotropic velocity, e.g., $U$, deviations from exterior values are radiated out at the speed of the external gravity waves:

$$U=U^{ext}-\sqrt{\frac{\mathrm{g}}{h}}\left(\zeta -{\zeta }^{ext}\right)$$

(10)

where ${\zeta }^{ext}$ and $U^{ext}$ are exterior values given from the model forcing along the boundaries (see below).

Furthermore, a radiation boundary condition was used for the 3D baroclinic velocities [53,54]:

$$\frac{\partial \phi }{\partial t}=-\left({\phi }_{\xi }\frac{\partial \phi }{\partial \xi }+{\phi }_{\eta }\frac{\partial \phi }{\partial \eta }\right)$$

(11)

where $\phi =\left\{u,v\right\}$, $\xi$ and $\eta$ are the curvilinear coordinates and

$$\left({\phi }_{\xi },{\phi }_{\eta }\right)=\frac{-\frac{\partial \phi }{\partial t}}{{\left(\frac{\partial \phi }{\partial \xi }\right)}^2+{\left(\frac{\partial \phi }{\partial \eta }\right)}^2}·\left(\frac{\partial \phi }{\partial \xi },\frac{\partial \phi }{\partial \eta }\right)$$

(12)

The barotropic timesteps were set to 5 min and 100 s for the NS and BCS models, respectively. Each barotropic step included 30 3D baroclinic timesteps without numerical instability problems occurring.

3.4. Model Forcing

The numerical model was forced with the combined action of astronomical tides and meteorological conditions. The open boundaries were forced using 20 tidal constituents, fifteen constituents were extracted from the TPXO9-atlas database [55]. The particular regional tidal solution used here as forcing factor relies on the European Shelf solution with a resolution of 1/30° [56]. The used dataset included 9 primary (M2, S2, N2, K2, K1, O1, P1, Q1, S1), 2 long period (Mf, Mm) and 4 non-linear (M4, MS4, MN4, 2N2) harmonic constituents. Five additional constituents (Ssa, L2, NU2, MU2, M6) were extracted from FES2014 with a resolution of 1/16° [57].

Atmospheric forcing in the present study used hindcast data obtained from the ECMWF (European Centre for Medium-Range Weather Forecasts) ERA5 Reanalysis data for the considered period [58]. The resultant ERA5 (ECMWF Reanalysis, fifth major global re-analysis) dataset is a combination of modelled and observed values, calculated using incremental 4D-Var data assimilation techniques via the ECMWF Integrated Forecast System (IFS), cycle 41r2 [58,59]. The ERA5-reanalysis database delivers published data on an hourly temporal scale and on a 0.25° × 0.25° spatial scale. The hydrodynamic model was run for the entire year 2011 to allow for full seasonal effects.

4. Model Validation

4.1. Observational Data

The numerical model was validated against observational sea level data measured during January 2011 to December 2011. Measurements from tidal gauges at 46 locations were used. The measuring points are shown in Figure 4. Eight of the tidal gauges lie within the Belgian coastal zone, which is of particular interest for the present study (see Figure 5).

4.2. Model Results for Astronomical Tides

The model is validated at the 46 measuring points using four statistical parameters: root-mean square error (RMSE), mean absolute error (MAE), correlation coefficient, and the skill index introduced by [60] and also used by [37] for ROMS efficiency evaluation:

$$S.I.=1-\frac{{\displaystyle \sum }{\left|X_{model}-X_{obs}\right|}^2}{{\displaystyle \sum }{\left(\left|X_{model}-\overline{X_{obs}}\right|+\left|X_{obs}-\overline{X_{obs}}\right|\right)}^2}$$

(13)

where $X$ is the variable being compared with a time mean $\overline{X}$. Perfect agreement between model results and observations will yield a skill of one and complete disagreement yields a skill of zero. The root mean square error is calculated from:

$$RMSE=\sqrt{\frac{{\displaystyle \sum }{\left(X_{model}-X_{obs}\right)}^2}{n}}$$

(14)

where $n$ is the number of observations.

Since most tidal models are forced by tidal constituents, it is a usual practice to validate model results against measured data in terms of the tidal constituents. Thus, in order to compare the model behavior with respect to the astronomical tides, a harmonic analysis was performed on both the measured and computed timeseries of the free surface elevation using T_Tide software [61]. Figure 6 depicts comparisons between the modelled and observed timeseries of the sea level due to astronomical tides for three of the Belgian tidal stations. The model results visually appear to be fairly well correlated with measurements.

As deducted from both literature and measurements throughout the Belgian Continental Shelf the semi-diurnal constituents M2, S2, N2 and K2 are dominant, with amplitudes higher than 15 cm [62]. There are also diurnal constituents, mainly O1 and K1, with significant energy, as well as higher order nonlinear constituents, such as M4 and MS4.

Table 1. Amplitude (m) of the three major tidal constituents in the Belgian Continental Shelf.

Stations	Observations			BCS Model			Deviations
	M2	S2	N2	M2	S2	N2	M2	S2	N2
Newport	1.92	0.57	0.33	1.94	0.57	0.34	0.02	0.00	0.01
Westhinder	1.69	0.50	0.29	1.68	0.49	0.29	−0.01	−0.01	0.00
Zeebrugge	1.66	0.47	0.28	1.69	0.49	0.29	0.03	0.02	0.01
Ostend	1.80	0.53	0.31	1.83	0.53	0.32	0.03	0.00	0.01
SW	1.58	0.43	0.29	1.65	0.47	0.28	0.07	0.04	−0.01
BVH	1.62	0.46	0.28	1.67	0.48	0.28	0.05	0.02	0.00
A2	1.65	0.47	0.28	1.69	0.49	0.29	0.04	0.02	0.01
Wandelaar	1.61	0.46	0.27	1.66	0.48	0.28	0.05	0.02	0.01

and

Table 2. Phase (degs) of the three major tidal constituents in the Belgian Continental Shelf.

Stations	Observations			BCS Model			Deviations
	M2	S2	N2	M2	S2	N2	M2	S2	N2
Newport	0.7	54.6	338.0	1.3	55.3	338.3	0.60	0.70	0.30
Westhinder	353.8	47.3	329.9	358.1	51.8	334.4	4.30	4.50	4.50
Zeebrugge	14.0	68.7	350.9	14.7	69.9	351.5	0.70	1.20	0.60
Ostend	0.5	54.4	337.7	6.8	61.3	343.7	6.30	6.90	6.00
SW	19.2	74.1	359.3	17.6	73.1	354.4	−1.60	−1.00	−4.90
BVH	14.8	69.8	352.0	15.1	70.4	351.9	0.30	0.60	−0.10
A2	11.5	66.1	348.6	13.1	68.1	349.8	1.60	2.00	1.20
Wandelaar	10.0	64.4	346.8	11.8	66.7	348.5	1.80	2.30	1.70

summarize the observed and modelled major tidal constituents for the eight Belgian tidal gauges used herein. For each station and constituent (M2, S2 and N2) the deviations between modelled and observed values of amplitude (in m) and phase (in degs) are also presented (modelled minus observed values).

Figure 7 depicts detailed modelled vs. measured data for the amplitude and phase of the 20 major tidal constituents at three Belgian tidal stations. The comparison between the observed and modelled data reveal a very good response from the ROMS model throughout the Belgian Continental Shelf. M2 is, by far, the most dominant constituent, having amplitudes that range from values slightly lower than 2.00 m to the west to approximately 1.60 m towards the northeast, near Scheldt estuary. The M2 tidal wave enters the Belgian coastal zone both from the north of the North Sea and the English Channel from the west. The model predicts this behavior with an average absolute bias of 4 cm and a tendency to slightly overpredict the amplitudes, with the exception of Westhinder station. The maximum of discrepancy of the amplitude is 7 cm and is observed at Scheur Wielingen. The phases of the eight stations are also well-predicted with a low average absolute error of 2.1°, which implies time lags of less than 5 minutes for the computed M2 tidal signal.

The second dominant constituent is the principal solar semi-diurnal S2 that mainly penetrates through the English Channel. The amplitudes vary from 0.40 m near the borders with Netherlands to 0.60 m near France. The model’s response is very good with respect to S2 resulting to an average absolute bias of less than 2 cm, while the phase deviations vary from 0.6° to 6.9°. N2 involves the ellipticity of moon’s orbit and its amplitude shows a restricted variation in Belgian waters, having values of approximately 0.30 m [62]. The model describes N2 amplitude with a significantly low error, while the phases deviate by a maximum 6.0°. As a comparison, [25] using TELEMAC 2D computed the amplitudes of M2, S2 and N2 at Wandelaar station with errors of approximately 9 cm, 3 cm and 7 cm, respectively whereas the ROMS model computes these constituents with errors of 5 cm, 2 cm and 1 cm, respectively.

Other significant semi-diurnal constituents in this region are the lunisolar K2, the smaller elliptic L2, the variational constituent MU2, the larger and the smaller lunar evectional constituents, NU2 and LDA2, respectively. These constituents have average over all stations amplitude that varies from 5 cm to 13 cm and the model can predict them with an average absolute error that does not exceed 1 cm. The agreement between modelled and measured data for the phases is also good, with the exception of the two last constituents at Scheur Wielingen station. The relevant absolute error is 16.6° for NU2 and 45.7° for LDA2, respectively.

Among the diurnal constituents (K1, O1, S1, P1, Q1) the most significant amplitudes are the ones of O1 and K1, with average measured values 0.10 m and 0.06 m, respectively. They are both accurately predicted by the model, with an accuracy of 1 cm. The corresponding maximum absolute errors for the phases are 10.9° for O1 at Westhinder station and 47.3° for K1 at Newport station, respectively. Furthermore, three long-period tidal constituents have been included in the analysis. The most important is the solar semiannual constituent (Ssa) with an almost uniform amplitude of 0.06 m throughout the Belgian coastal zone. The lunar monthly constituent (Mm) has average amplitude of 0.04 m, while the lunar fortnightly (Mf) is negligible (~1 cm). All three constituents are accurately predicted by the model near Belgium with an error of less than 1 cm. With regard to the phases of the lunar long-period constituents, the minimum errors are observed at Westhinder (0.5° for Mm) and Newport (14.7° for Mf), respectively, while the maximum at Scheur Wielingen (30.5° for Mm and 35.0° for Mf, respectively). On the other hand, the phase of Ssa is computed with a minimum absolute error of 1.0° at A2 station and a maximum of 13.3° at Newport.

As the tidal waves propagate into the shallow water of the Belgian Continental Shelf higher-order nonlinear components are generated from the interactions of primary constituents. The ones with greatest amplitude in this region of the North Sea are M4, M6, MS4, MN4, and 2N2. The largest shallow water constituent in the northwest European Shelf is M4. This constituent is generated by self-interactions of M2 constituents in shallow water. Mathematically, the origin of M4 is due to time-varying depth and the advection terms and is known for causing double high and low waters along parts of the Dutch and British coasts in the English Channel [63]. When double high water occurs, only a few hours remain for the low water. Due to this short duration of low water very strong currents are formed. Throughout the Belgian coastal zone the measured amplitude of M4 varies from 0.09 m to 0.14 m, with the largest values observed near Newport. The model computes these values very accurately, with absolute errors of less than 1 cm. As far as the phases are concerned, the absolute error has an average value of 8.9° and a maximum of 17.7° at Ostend. In contrary to the advective origin of M4, the shallow water overtide M6 results from M2 self-interactions due to bottom friction [64]. This sixth-diurnal constituent is particularly important in the region of Le Havre in the English Channel, having measured amplitudes that may exceed 15 cm. Within the Belgian territory its average amplitude is 8 cm, and it may even reach 10 cm near Scheur Wielingen. The model underestimates these values, giving an average absolute error of 0.03 m in Belgium. Increased discrepancies are also observed with respect to the phases, which are computed with an average absolute error of 36.2°. The deviations may even be as high as 99.4° at Newport station. On the other hand, the amplitude of the less significant shallow water constituents MS4, MN4 and 2N2 is very accurately predicted by the model, with errors of less than 1 cm. Nevertheless, deviations between the modelled and observed phases are again present. The highest of all are for 2N2 constituent, for which the average absolute error of Belgian stations is 70.8°.

Since the model is forced with a limited number of tidal constituents, the above constituent-by-constituent analysis offers a better evaluation of the model performance than the direct metrics of the entire timeseries of water level, e.g., MAE, RMSE. In this context, the average absolute error of Belgian stations with respect to amplitude is 2.3% for M2, 3.7% for S2 and 2.4% for N2, respectively. The higher values are observed for M6 (43.2%) and MF (25.2%) constituents. The aforementioned statistical parameters for the timeseries of water level due to astronomical tides were also calculated and they are summarized in

Table 3. Statistical parameters and error with regard to astronomical tides and atmospheric forcing at eight Belgian stations and various other stations in the North Sea.

Tidal Stations	Astronomical Tides				Atmospheric Forcing
	MAE (m)	RMSE (m)	Cor	S.I.	RMSE (m)
22. Newport	0.081	0.100	0.998	0.999	0.114
23. Westhinder	0.104	0.121	0.996	0.998	0.110
24. Zeebrugge	0.073	0.091	0.998	0.999	0.114
25. Ostend	0.156	0.179	0.992	0.996	0.118
26. Scheur Wielingen	0.089	0.116	0.997	0.998	0.118
27. Bol Van Heist	0.073	0.090	0.998	0.999	0.125
28. A2	0.080	0.100	0.997	0.998	0.121
29. Wandelaar	0.082	0.102	0.997	0.998	0.126
46. Stavanger	0.026	0.033	0.969	0.984	0.081
44. Hanstholm	0.031	0.037	0.947	0.954	0.051
42. Hörnum	0.112	0.132	0.993	0.989	0.075
35. Platform A12	0.081	0.098	0.894	0.944	0.133
20. Calais	0.210	0.259	0.991	1.000	0.196
11. Bournemouth	0.181	0.221	0.829	0.908	0.153
4. Whitby	0.128	0.151	0.959	0.996	0.203
2. Aberdeen	0.129	0.161	0.989	0.994	0.181

. The high values of correlation and skill index (S.I.), near unity, reveal a very good response of the model for the Belgian coastal zone. The best results are at Bol Van Heist and Zeebrugge with mean absolute error of less than 7.3 cm. The average value of MAE for the Belgian stations is 9.2 cm. The corresponding average RMSE is 11.2 cm, with the lowest values resulting at Bol Van Heist and Zeebrugge stations, i.e., approximately 9.0 cm. On the other hand, the highest errors are observed in Ostend. The skill index is very close to unity for all Belgian stations, from 0.996 to 0.999, revealing a very good degree of agreement between observations and model results.

The results for eight additional stations around the North Sea are also depicted in Table 3. The model simulations are deemed to be very good at the northeast, i.e., Denmark and Norway, with RMSE less than 4 cm. The accuracy was also found good at the German Bight and near the Dogger Banks with RMSE less than 0.14 m. Indicatively, the RMSE at Hörnum station (Sylt island, northwest Germany) and the Dutch Platform A12 (central North Sea) was 0.132 m and 0.098 m, respectively. The accuracy is also good along the east coast of the United Kingdom, with the RMSE at Aberdeen and Whitby being 16.1 cm and 15.1 cm, respectively. On the other hand, the model’s behavior in the English Channel, especially along the French coasts in the north of Normandy, is poorer than in other parts of the domain and the RMSE may exceed 25 cm. Perhaps, a spatially variable bed friction coefficient could result to reduced values of the RMSE. This is an issue that can be further investigated. However, even in this area the skill index is always higher than 0.90, which reveals an overall adequate response of the model [37].

Figure 8 shows the co-tidal map of the M2 constituent as computed both from the model and the observations. The location of the two amphidromic points, in the Southern Bight and the Dogger Banks, respectively, is predicted very accurately by the model. The third one, southwest of the Norwegian coast, is slightly shifted to the southeast. The same co-tidal map was also extracted from the TPXO9-atlas database, and the model results for both amplitudes and phases seem to be in good agreement with this benchmark database. As described above, Kelvin waves penetrate from the north and they are strongly dissipated by bottom friction in the shallow southern coastal waters [29]. As a result, increased amplitudes are observed towards the eastern coast of the United Kingdom and significantly smaller amplitudes off the coasts of Denmark and Norway. These results can also be extracted from Figure 8 for both the model and the TPXO9 database.

4.3. Model Subtidal Results

The principal driver of subtidal circulation in the North Sea is the wind forcing, and secondarily the local gradients of air pressure [33]. The model internally computes the ocean-atmosphere boundary layer based on [66]. This module relies on a number of analytical formulations which take into account several parameters, e.g., surface air temperature, surface air relative humidity, cloud fraction, evaporation, precipitation, short- and long-wave radiation fluxes, etc. The resulting boundary conditions refer to the surface air pressure provided as external forcing and the surface stress due to wind (and rainfall) given by:

$${\mathit{\tau }}_{\mathit{s}}={\rho }_{air}·C_d·\sqrt{{\left|{\mathit{U}}_{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}\right|}^2+w_{gus}{}^2}·{\mathit{U}}_{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}+0.85·{\phi }_{rain}·\left|{\mathit{U}}_{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}\right|·\left(\frac{U_{wind}}{\left|U_{wind}\right|},\frac{V_{wind}}{\left|V_{wind}\right|}\right)$$

(15)

where ${\mathit{U}}_{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}=\left(U_{wind},V_{wind}\right)$ is the wind velocity 10 m above the MSL. The contribution of gustiness to wind speed is expressed through the term $w_{gus}$, while ${\phi }_{rain}$ (=0 here) is the rainfall rate, ${\rho }_{air}$ is the air density, and $C_d$ is a parameter based on Monin-Obukhov similarity theory [66]. The wind gustiness is difficult to be predicted and, in general, it is more important in coastal areas. However, its effect has been found to be not significant in the North Sea [67].

In order to validate the model’s efficiency in reproducing the hydrodynamics induced by the varying atmospheric conditions two simulations were conducted, one with and one without the inclusion of meteorological forcing, i.e., wind and air pressure fields. By subtracting the latter from the former, the aforementioned atmospheric driver can be isolated and studied. As far as the observations are concerned, a Lanczos-cosine low-pass filter was applied on the measured timeseries in order to calculate the effect of atmosphere on the circulation [68]. The resulting timeseries from the model results and the observations were compared with each other. This comparison for three of the Belgian stations is shown in Figure 9.

The resulting errors from these comparisons are summarized in Table 3. The RMSE error is calculated from the modelled timeseries and the observed timeseries, after applying the aforementioned low-pass filter to the latter. The former results from the difference of the model output when the atmospheric forcing is excluded from the corresponding values which include the effect of atmospheric forcing. The average RMSE in the Belgian coastal zone is 11.8 cm, with the minimum value at Westhinder (11.0 cm) and the maximum error at Wandelaar stations (12.6 cm). The correlation between the modelled and the observed water levels is shown in Figure 10. The correlation coefficient varies from $r=0.67$ at Westhinder to $r=0.76$ near the estuary of Scheldt River. The RMSE for some additional stations around the North Sea is also presented in Table 3. Its values vary between approximately 5–13 cm in the central and northeast domain and they increase towards the eastern coast of the United Kingdom, where values of 20 cm are calculated. The RMSE in the English Channel lies within, approximately, the range 15 to 20 cm. In general, the results are comparable, and in some cases better than other studies. [33] presented results of 5 different models simulating storm events in the North Sea with a duration of a few days. They computed bias at Ostend station that may be as high as 35 cm and standard deviations varying from 10 to 42 cm. In the present study, the corresponding RMSE over the entire 2011 year was 11.8 cm. In Aberdeen and Stavanger, the standard deviation for the various models studied by [33] varied in the range 6–19 cm and 3–10 cm, respectively. Herein, the corresponding RMSE are 18.1 cm and 8.1 cm, respectively.

4.4. Thermohaline Variation and Atmospheric Surface Fluxes

It is stated above that the southern part of the North Sea, and particularly the Belgian waters, constitute a well-mixed environment. That is why in the described, so far, configuration the water temperature and salinity were assumed uniform and time invariant. However, this does not apply for the entire North Sea. In order to check the model’s efficiency in reproducing the thermohaline variations in the North Sea, the model was also applied including the variation of the water temperature and salinity. The initial conditions for these two variables were extracted from the ECMWF ERA5-reanalysis database and a nonlinear equation of state, Equation (6), was employed.

In order to further increase the realism of the model configuration, the bulk parameterization used in the previous sections [48] is here replaced by the full atmospheric forcing based on meteorological data provided by ECMWF ERA5-reanalysis database [58,59]. In this context, the values of the surface wind stress, ${\mathit{\tau }}_{\mathit{s}}=\left({\tau }_{s,x},{\tau }_{s,y}\right)$, were directly extracted from ECMWF database, and they were not computed internally by the model through Equation (15). However, it should be noted that the calculation of the wind stresses in ERA5 relies on an ocean state that might be different from the ocean state in ROMS. This fact may result to mismatches between the modelled and the imposed air-sea interface, respectively. Additional variables that were obtained from the global database and included in the atmospheric forcing of the model were: (a) the surface net solar (shortwave) radiation flux, (b) the surface net thermal (longwave) radiation flux, (c) the surface sensible heat flux, (d) the surface latent heat flux, (e) the rain precipitation rate, and (f) the evaporation rate. The required by ROMS net surface heat flux is given from the sum of the first four variables.

A comparison between the modelled and the observed annual average Sea Surface Temperature (SST) for 2011 is depicted in Figure 11. The observations were taken from ECMWF ERA5—reanalysis database. From Figure 11 it can be claimed that the model can qualitatively predict the variation of the annual average SST over the entire North Sea. The good agreement between observations and model results is confirmed by a relatively low RMSE of 0.68 °C over the entire domain. The maximum difference between the observations and the model results is 2.26 °C. The figure also confirms the conclusion that within the Belgian Continental Shelf the water is well-mixed and the maximum deviation between the model and the observations is 1.22 °C.

Figure 12 shows the vertical distribution of the annual and spatial (over the entire North Sea) average potential temperature, as this is both extracted from GLORYS12V1 reanalysis database and computed by the presented model. The agreement is fairly good for the first 40 m of depth, where the deviations are less than 0.5 °C. In the zone between 50 and 70 m of depth, the model overpredicts the water temperature, with the deviations ranging between (0.55–0.70) °C. Deeper than this zone, up to approximately 125 m of depth, the agreement is again very good, with differences less than 0.2 °C. In deeper layers the model constantly underpredicts the temperature, with the maximum deviations being of 1.25 °C at approximately a depth of 200 m. Throughout the considered domain, these great depths (>200 m) are only met along the Norwegian Trench. From both curves, it is observed that the average temperature at the water surface is approximately 10.3 °C and drops to 6.8 °C over the first 100 m of depth. Up to the depth of 500 m the average temperature has further dropped to 4 °C.

Forcing the model with the realistic atmospheric data from ECMWF reanalysis database has also resulted to an improvement of its behavior with respect to meteorological conditions.

Table 4. RMSE and correlation coefficient with regard to atmospheric forcing at eight Belgian stations and various other stations in the North Sea when applying: (a) the bulk formulation and (b) data from ECMWF reanalysis for meteorological variables.

Tidal Stations	Atmospheric Forcing—Bulk Formulation [48]		Atmospheric Forcing—ECMWF
	RMSE (m)	Cor	RMSE (m)	Cor
22. Newport	0.114	0.707	0.085	0.762
23. Westhinder	0.110	0.668	0.085	0.696
24. Zeebrugge	0.114	0.754	0.106	0.812
25. Ostend	0.118	0.730	0.084	0.784
26. Scheur Wielingen	0.118	0.587	0.094	0.653
27. Bol Van Heist	0.125	0.756	0.104	0.813
28. A2	0.121	0.729	0.094	0.790
29. Wandelaar	0.126	0.740	0.110	0.790
46. Stavanger	0.081	0.060	0.031	0.122
44. Hanstholm	0.051	0.032	0.039	0.045
42. Hörnum	0.075	0.319	0.029	0.885
35. Platform A12	0.133	0.086	0.032	0.635
20. Calais	0.196	0.173	0.195	0.173
11. Bournemouth	0.153	0.172	0.041	0.409
4. Whitby	0.203	0.220	0.088	0.451
2. Aberdeen	0.181	0.015	0.180	0.099

summarizes the RMSE and correlation coefficient using these forcing data, in comparison with the corresponding values when applying the surface bulk flux formulation. The average reduction of the RMSE within the BCS is 2.3 cm, which implies an improvement of about 20%. The maximum improvement is observed at Ostend, with a reduction of the RMSE of 3.4 cm.

The RMSE and correlation for some additional stations around the North Sea are also presented in Table 4. In the central and eastern subdomain the RMSE varies between 3–4 cm, which consists a reduction of 25–75% compared to the case of using the bulk flux formulation. An improved behavior is also noticed towards the eastern coast of the United Kingdom, e.g., a reduction of RMSE of 9 cm at Whitby. Towards the northeast coast of U.K. this improvement is reduced. For example, at Aberdeen the RMSE is almost identical, but the correlation is improved (0.099 instead of 0.015 using the bulk flux formulation). Within the English Channel, the model response is significantly improved along the coast of the United Kingdom, e.g., at Bournemouth a reduction of the RMSE of 73%. However, along the French coast the higher values of error remain unchanged.

A comparison between the observed timeseries of the water level due to purely atmospheric forcing and the model results using the two formulations described above, respectively, is shown for three Belgian stations in Figure 13. The improvement of the model behavior is also confirmed with respect to the correlation coefficient, as shown both in Table 4 and Figure 14. The average correlation coefficient within the Belgian coastal zone is now $r=0.763$, as compared to $r=0.709$ when using the bulk flux formulation. Throughout the entire North Sea, the improvement with respect to correlation can be as high as 0.55, at the central and northeast domain (stations Hörnum and Platform A12 in Table 4).

5. Discussion and Conclusions

A hydrodynamic model for the NW European Shelf was developed based on ROMS and its application to the North Sea and the Belgian Continental Shelf was discussed. The flexibility of ROMS, which allows the user to switch on/off a particular physical process, is a major advantage for this model. However, in large-scale applications it is always a challenge to balance an accepted accuracy level with the simplicity resulting from selecting global physics modules applicable over the entire domain. For example, the southern part of North Sea is a well-mixed environment, which is not the case for the northern parts, especially in summer [29]. As a result, the thermohaline stratification was neglected in the model as the area of particular interest is the Southern Bight in addition to the sufficiently strong tidal currents that retain the fully mixed water column even when deeper, offshore waters are stratified [69].

In the Belgian Continental Shelf the model succeeds in describing the amplitude of the dominant harmonic constituents M2, S2 and N2, with absolute errors of 2.3%, 3.7% and 2.4%, respectively. The accuracy of the predicted diurnal and long-period constituents is also fairly good. Thus, the major part of the tidal energy within the Belgian coastal zone is accurately described. In coastal areas some nonlinear shallow-water constituents may have increased amplitude, as they carry non-negligible energy portions due to the shallowness of the North Sea in general. Particularly in the Southern Bight, M4 and M6 are the major higher-order constituents. M4 is by far the most important and the model computes very accurately its amplitude with an average absolute error of 1 cm, while the phase has an average absolute error of 22.7° in the Belgian waters. On the other hand, significant errors, that may reach 50–60%, are observed with respect to M6. Accurate description of these higher-order constituents is important for operational reasons in the coastal zone, as they are related to tidal asymmetry and correct prediction of high and low waters. The aforementioned relative errors correspond to an average absolute error for M6 of 3 cm which is not large enough, alone, to significantly alter the prediction of double high and low waters.

The overall tidal pattern is correctly simulated by the model. The three amphidromic points of M2 constituent are well predicted. This is also true for the Kelvin tidal wave along the eastern coastline of the United Kingdom. In the Norwegian Trench and western Jutland the model is very accurate, with RMSE due to astronomical tides of less than 4 cm. The model accuracy is also good in the central North Sea and the German Bight having RMSE of approximately 10 to 15 cm. On the other hand, significant inaccuracies are observed in the English Channel with errors that may locally exceed 25 cm. However, these errors are reduced towards the Belgian and Dutch coast and the Southern Bight, as the result of the very shallow bathymetry and the combined action of tidal waves coming from both west and north.

The model’s efficiency in reproducing the atmospheric induced hydrodynamics was also tested. The model was run with and without atmospheric forcing and the difference between the two simulations gives the effect of meteorological conditions on the model results. The correlation between measurements and model results within the Belgian coastal zone for a model that includes only the atmospheric forcing is weaker than in the case of where only the tides (model neglecting atmospheric forcing) is considered and varies between $r=0.666$ and $0.755$, as opposed to $r=0.992$ and $0.998$, while the average RMSE is 11.8 cm, as opposed to RMSE = 11.2 cm for the model without atmospheric forcing included. Similarly to the tides, the model behaves quite well to the northeast of the domain and the central North Sea, with RMSEs between 5 and 13 cm. Along the eastern British coast the inaccuracy of the model is increased both compared to the interior of the domain and the corresponding tidal-induced water level variation. On the other hand, the error in the English Channel is computed somewhat lower, between 15 and 20 cm.

The reasoning for the variation of the model skill and, particularly, the lower accuracy in the English Channel, rely on a combination of factors. At first, it is the relatively steep seabed there and the rather poor resolution of 5 km × 5 km to resolve these slopes. Secondly, as far as the open boundary conditions are concerned, only the tidal forcing (harmonic constituents) from TPXO9 was imposed, and not full timeseries of the involved variables downscaled from a global (or of a larger scale) model. Thirdly, the spatial and temporal variations in bottom roughness were neglected as a constant friction coefficient was used. Fourthly, the observations were collected from different sources for the various countries and, thus, their accuracy may also vary. Finally, despite using a reliable database for the atmosphering forcing, the atmospheric models always contain a degree of uncertainty, which in combination to the level of accuracy of the modelled surface boundary layer, may have an effect on the accuracy of the ocean model.

The model’s behavior with respect to thermohaline variations was also investigated. It was concluded that it can accurately reproduce the variations of the sea surface temperature throughout the entire North Sea, and with even higher accuracy within the Belgian coastal zone. The vertical variation of water temperature was also predicted with good accuracy, particularly in the upper 125 m from the sea surface.

The importance of using realistic forcing data for the most significant atmospheric variables (in addition to the wind velocities and air pressure) was also investigated. By forcing the model with data from the ECMWF global database with respect to wind stresses and surface fluxes, both the RMSE and the correlation were significantly improved, as compared to using a surface layer bulk flux formulation. However, in order to draw a general conclusion about the superiority of using either the ERA5 full forcing or the bulk flux formulation, it is necessary that all the required variables in the latter (humidity, rain, air temperature, radiation) are varying in space and time.

In conclusion, the research questions posed in the introduction have been largely addressed. A new hydrodynamic model has been independently developed with particular focus on maximizing its accuracy in the Belgian Continental Shelf. This included a more accurate set of bathymetric data, a finer grid resolution of the model and the calibration with respect to a number of parameters, i.e., bottom friction coefficient and parameters of the GLS scheme. The model has been validated both for astronomical tides and atmospheric forcing and showed very good response in the Belgian waters. It can also be used for validating other numerical models. Based on a limited number of available results in the literature, the present model shows similar or better accuracy with respect to principal or high-order tidal constituents and atmospherically induced water level variations. In addition, it can be concluded that the wide range of options and modules included in ROMS offer the opportunity to use the model for more focused research on the accuracy of the used input and simulated processes. Finally, the most important added value of this work is that a validated hydrodynamic ROMS model forms a basis for further study using a research type model for various scientific or operational applications.