A Regional Operational Model for the North East Atlantic: Model Configuration and Validation

Abstract

An operational model for an area of the northeast Atlantic that encompasses all of Ireland’s territorial waters has been developed. The model is an implementation of the Regional Ocean Modelling System (ROMS) and uses operationally available atmospheric and boundary forcing, and a global tide solution for tidal forcing. River forcing is provided by climatological daily discharge rates for 29 rivers across Ireland, west Britain, and west France. It is run in an operational framework to produce 7-day hindcasts once a week, and daily 3-day forecasts which are published in a number of formats. We evaluated the model skill by comparing with measured data and calculating statistics such as mean error, root mean square error (RMSE), and correlation coefficient. The observations consist of satellite Sea Surface Temperature (SST), total surface velocity fields from satellite, water level time series from around the Irish coast, and temperature and salinity data from Array for Real-Time Geostrophic Oceanography (ARGO) and Conductivity Temperature Depth (CTD) profiles. The validation period is from 1 January 2016 until 31 December 2019. The correlation coefficient between the model and satellite SST is 0.97 and recorded in March and April 2018. The model error is about 5% of the total M₂ amplitude in the Celtic Sea recorded at Dunmore East tide gauge station. The maximum RMSE between the model and the CTD temperature profiles is 0.8 °C while it is 0.17 PSU for salinity. The model correctly defines the shelf water masses around Ireland. In 2019 the Irish Coastal Current (ICC) was very strong and well defined along most of the western Irish coast. The model results have well reproduced the ICC front for the whole simulation period.

1. Introduction

Ireland is located in a hydrographically interesting part of the northeast Atlantic Ocean where the regimes of the deep ocean and the continental shelf interact. The North Atlantic Current (NAC) divides southwest of the Rockall Bank and produces several branches flowing roughly to the northeast or southeast [1]. The Rockall Trough is one pathway by which heat is transported into the Nordic Seas [2,3] (see Figure 1). The continental shelf edge has a poleward current associated with it, which is seasonally variable [4,5,6] and is the interface for the ocean-shelf exchange of nutrients and carbon [7]. On the Irish shelf, the Irish Shelf Front (ISF) separates fresher Irish coastal waters from the oceanic Eastern North Atlantic Water (ENAW) at ca. 11° W [8]. The Shannon River on the west coast of Ireland is the largest river in the region with a mean annual discharge of 186 m³ s⁻¹ [9], but there are also significant freshwater discharges from rivers on the south coast of Ireland and west coast of Britain as well as from further afield such as the Loire in France. This induces vertical density stratification on the shelf while thermal stratification is an important process during the summer months. In contrast, mixing due to tide and wind creates areas with well-mixed water, such as the Irish Sea. Fronts occur between these stratified and well-mixed areas, such as the Celtic Sea front, which is located at the southern entrance to the Irish Sea [10]. A seasonal density-driven current, the Irish Coastal Current (ICC) has been described by [10,11]. The ICC is linking Land’s End in southwest England to the northwest coast of Ireland and this may provide an important pathway for harmful phytoplankton, which has significant consequences for the aquaculture industry on the west coast of Ireland and onwards to the coast of west Scotland. The Scottish Coastal Current flows northwards from the Irish Sea along the west coast of Scotland [12].

Operational modelling of physical ocean fields has developed rapidly over the last decade in response to the needs of the marine community. The availability of hindcasts and forecasts of the hydrodynamic properties for specific regions is of great benefit to decision makers in areas such as marine environmental management, and the aquaculture and fisheries industries. Operational models that forecast the physical state of the oceans exist for European oceans (e.g., [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]). The first operational models developed in Ireland and described by [15,16].

Ref. [24] described an operational biogeochemical model for the North East Atlantic Ocean region. The model domain covered a significant portion of the North-West European continental shelf with a horizontal and vertical resolution of ≈4 km and 40 sigma levels, respectively. It was one-way nested within the high resolution (1/12°) CMEMS Ocean (http://www.CMEMS-ocean.fr) PSY2V4R2 operational model of the North Atlantic whereby daily values for potential temperature, sea surface height and velocity are linearly interpolated at the open ocean boundaries. A simple nitrogen-based NPZD model developed described in [28] was dynamically coupled to ROMS. The model has successfully represented the intra-annual variability of surface chlorophyll and nitrate concentrations at monthly time scales.

There have already been several studies simulating the North Atlantic region including [27,29,30,31,32,33,34,35,36]. Overall, properly validated hydrodynamic models can provide information on the present and near-future states of ocean systems, which is of value to a wide range of users from both the private sector and the general public.

Objectives of the Study

This study presents the set-up and results from the validation of a high-resolution numerical operational model (hereafter called NEA_ROMS) developed at Irish Marine Institute. Validation of NEA_ROMS was carried out for the time period from January 2016 until December 2019 using observational data from various sources. The model has been compared with satellite SST, tide gauges, Argo, and CTD temperature and salinity profiles. In addition, the Geostrophic and Ekman Current Observatory GEKCO surface data was used to validate model velocity fields. Additionally, this study will focus on the Irish coastal waters and the authors will examine the representation of the Irish coastal current (ICC) in the model. Section 2 describes the model implementation, and nesting procedures, Section 3 presents the validation against observational data, Section 4 describes the model results related to the ICC, and Section 5 provides the discussion and conclusions.

2. Model Design and Implementation

The model used in this study is an implementation of the Regional Ocean-Modelling System (ROMS). ROMS has been proven to demonstrate substantial skill in forecasting [37] and has been used in a number of successful operational systems [38,39,40]. ROMS is a primitive-equation, free-surface, hydrostatic model as described in [41]. In the horizontal, it uses orthogonal curvilinear coordinates on an Arakawa “C” grid while utilizing a terrain-following (sigma) coordinate in the vertical. The momentum equations are solved using a split-explicit time-stepping scheme whereby a whole number of barotropic time steps are carried out within each baroclinic time step to solve the barotropic momentum equations. The prognostic variables of the model are surface elevation, potential temperature, salinity, and velocity. All the model equations are written in rectangular coordinates and contain spatially and temporally varying horizontal eddy viscosity and diffusion coefficients. Vertical diffusion is calculated using the K-profile parameterization (KPP) [42] and a modified Jacobian method [43] is used to calculate horizontal pressure-gradient forces. The associated parameters with KPP vertical turbulent closure scheme have been enhanced to tune the vertical profile of currents, temperature, and salinity by [44,45,46,47]. Recursive MPDATA is used for the advection of tracers [48]. Bottom stress is applied using the logarithmic “law of the wall” with a constant roughness length of 0.01 m.

The NEA_ROMS model domain covers a significant portion of the North-West European continental shelf and also the Porcupine and Rockall Banks and the Rockall Trough as shown in Figure 1. NEA_ROMS model has a horizontal resolution of 1.1 km and 40 vertical sigma levels with the highest concentration of levels at the surface and the bottom. The basin-scale daily average lateral boundary conditions are provided from the high resolution (1/12°) CMEMS_GLOBAL_ ANALYSIS_ FORECAST_ PHY _001_024 model [49]. The daily values are for potential temperature, salinity, sea surface height, and velocity components at the three open boundaries. The discharges from 29 rivers are included in the NEA_ROMS model. The rivers comprise the major rivers of Ireland, west Britain, and west France, as well as some more minor Irish rivers. The discharge rates are average daily values calculated from many years of historical discharge data for 29 rivers across Ireland, west Britain, and west France. Furthermore, the nested model temperature and salinity fields were nudged toward the GLOBAL_ ANALYSIS_ FORECAST_ PHY _001_024 operational model. Within the NEA_ROMS, several one-way nested shelf models have been implemented to downscale and forecast near the Irish coasts.

The model bathymetry utilizes data from a number of sources. Large-scale multi-beam mapping of Irish territorial waters has been carried out under the Irish National Seabed Survey (INSS) and the Integrated Mapping for the Sustainable Development of Ireland’s Marine Resource (INFOMAR) programs. These datasets along with an extensive single-beam archive maintained by the Marine Institute were used to develop the model bathymetry. A 1 × 1 km Celtic Sea bathymetry [50] complements the bathymetry in the Irish and Celtic Seas while GEBCO data is used for the bathymetry in non-Irish waters. The model minimum depth is set at 10 m near the coastline. In order to reduce the pressure gradient errors, we smoothed the model bathymetry according to the rx₀ factor of Beckman and Haidvogiel [51,52].

2.1. Lateral Boundary Conditions

The model has three open boundaries in the North, East, and West. The NEA-ROMS lateral boundary conditions are taken from the CMEMS GLOBAL_ANALYSIS_FORECAST_ PHY _001_024 at 50 vertical levels are ranging from 0 to 5500 m. The CMEMS model is version 3.1 of NEMO [49]. The model data is freely available via the Copernicus Marine Environment Monitoring Service (CMEMS) web site marine.copernicus.eu. The analysis is updated weekly while a 10-day forecast is updated daily. The model boundary parameters are temperature, salinity, total velocity components, surface elevation, and barotropic velocity components. The nesting was designed to ensure that the volume transport across the open boundary of the nested model matches the volume transport of the CMEMS global model.

For the barotropic velocity, we used the scheme proposed by [53]. At the outflow, we impose

$$\overline{{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}}=\left[\sqrt{\frac{\mathrm{g}}{{\mathrm{H}}_{\mathrm{nested}}}}\left({\mathsf{\eta }}_{\mathrm{nested}}-{\mathsf{\eta }}_{\mathrm{nesting}}\right)\right]-\overline{{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}}\frac{{\mathrm{H}}_{\mathrm{nesting}}+{\mathsf{\eta }}_{\mathrm{nesting}}}{{\mathrm{H}}_{\mathrm{nested}}+{\mathsf{\eta }}_{\mathrm{nested}}}$$

(1)

where $\overline{{\mathsf{\phi }}_{\mathrm{n}}}$=$\frac{1}{\mathrm{H}+\mathsf{\eta }}{\displaystyle \underset{\mathrm{H}}{\overset{\mathsf{\eta }}{\int }}{\mathsf{\phi }}_{\mathrm{n}}\mathrm{dz}}$ is the normal barotropic velocity at the open boundaries and g the gravity, while at the inflow will be $\overline{{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nesting}}}=\overline{{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}}$

For the baroclinic velocity, we impose the normal velocity according to the following formula:

$${\displaystyle \underset{{\mathrm{n}}_2}{\overset{{\mathrm{n}}_1}{\int }}{\displaystyle \underset{{\mathrm{H}}_{\mathrm{nested}}}{\overset{{\mathsf{\eta }}_{\mathrm{nested}}}{\int }}{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}}}\mathrm{dzdn}={\displaystyle \underset{{\mathrm{n}}_2}{\overset{{\mathrm{n}}_1}{\int }}{\displaystyle \underset{{\mathrm{H}}_{\mathrm{nesting}}}{\overset{{\mathsf{\eta }}_{\mathrm{nesting}}}{\int }}{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nesting}}}}\mathrm{dzdn}$$

(2)

where ${\mathsf{\phi }}_{\mathrm{n}}$ is the normal component to the three open boundaries in m/s. For the temperature and salinity, we used a third-order upstream scheme with implicit mixing for horizontal advection [54]. Whenever the normal velocity is directed outwards:

$$\frac{\mathsf{\partial }{\mathrm{T}}_{\mathrm{nested}}}{\mathsf{\partial }\mathrm{t}}+{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}\frac{\mathsf{\partial }{\mathrm{T}}_{\mathrm{nested}}}{\mathsf{\partial }\mathrm{y}}=0;\frac{\mathsf{\partial }{\mathrm{S}}_{\mathrm{nested}}}{\mathsf{\partial }\mathrm{t}}+{{\mathsf{\phi }}_{\mathrm{n}}}_{\mathrm{nested}}\frac{\mathsf{\partial }{\mathrm{S}}_{\mathrm{nested}}}{\mathsf{\partial }\mathrm{y}}=0$$

(3)

While on the inflow the temperature and salinity are prescribed to be identical to that interpolated to CMEMS global model values. Where n1, n2 are the extremes of the open boundary section, ηnesting, and Hnesting are the surface elevation and the bathymetry of the nested model, while ηnested and Hnested are the surface elevation and bathymetry of the CMEMS global model and ${\mathsf{\phi }}_{\mathrm{n}}$ is the normal velocity component to the open boundaries. The total tangential velocity at the open boundaries is set to be equal between the CMEMS global model and the nested model.

The nested model temperature and salinity fields were nudged toward the CMEMS global model values in the right-hand side (r.h.s) of the prognostic fields following [55]. Indicating temperature or salinity with ‘gamma’ we add a relaxation term such as:

$$\frac{\mathsf{\partial }\mathrm{\gamma }}{\mathsf{\partial }\mathrm{t}}=\mathrm{r}.\mathrm{h}.\mathrm{s} \frac{1}{\Gamma }\left(\mathrm{\gamma }-{\mathrm{\gamma }}^{\mathrm{CMEMS}}\right)$$

(4)

Here Γ is the relaxation time that varies smoothly from 20 model grid internal points to the open boundaries with values of 240 s to 100 days. This nudging prevents a substantial model drift from the CMEMS global model values.

Finally, the tidal forcing is prescribed at the model boundaries by applying elevations and barotropic velocities for 10 major tide constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁, M_f, M_m) which are obtained from the TPXO8 global inverse barotropic tide model [56]. Radiation conditions are specified at the boundaries for depth-averaged velocity after [53] and also for water level after [57]. Radiation conditions are prescribed for tracers and baroclinic momentum [55].

2.2. Surface Boundary Conditions

The atmospheric data for the computation of the surface forcing are taken from the global high resolution (0.125°) atmospheric model run by the European Centre for Medium Range Weather Forecasts (ECMWF) with 3-hourly temporal step. This has changed in 2019 to hourly fields with spatial resolution (0.1°). The atmospheric fields used are air temperature, relative humidity, wind velocity at 10 m above sea level, mean sea level pressure, cloud cover, total precipitation, and solar radiation (net surface and longwave radiation). In order to parametrize the air–sea interaction processes, the wind stress, heat fluxes, and evaporation rate are computed by means of interactive bulk formulae making use of atmospheric data and the model predicted sea surface temperature according to [58]. The NEA-ROMS heat fluxes are computed interactively according to [13] and described in [59] the heat flux boundary condition is given by

$${\mathsf{\rho }}_{\mathrm{o}}{{\mathrm{C}}_{\mathrm{P}}{\mathrm{K}}_{\mathrm{H}}\frac{\mathsf{\partial }\mathrm{T}}{\mathsf{\partial }\mathrm{z}}|}_{\mathrm{z}=\mathsf{\eta }}={\mathrm{Q}}_{\mathrm{T}} +\frac{\mathsf{\partial }\mathrm{Q}}{\mathsf{\partial }\mathrm{T}}{ (\mathrm{SST}}_{\mathrm{CMEMS}}{-\mathrm{SST}}_{\mathrm{m}})$$

(5)

where Q_T is the total net heat flux according to the formula:

Q_T = Q_S − Q_b − Q_e − Q_h

(6)

The last term in Equation (5) is the heat flux correction term. ∂Q/∂T is the heat flux correction. which is set at 80 W·m⁻²·°C⁻¹ for NEA-ROMS. SST_m is the model sea surface temperature and salinity and SST_CMEMS is the reference SST. Q_S, Q_e, Q_b, and Q_h are the short, latent, back, and sensible heat fluxes in W/m². (K_H) is the vertical mixing coefficients for tracers in (m²·s⁻¹). This heat flux correction, which was extracted from CMEMS, sea surface temperature, is to control the simulation from drifting away from the CMEMS values.

The momentum flux boundary condition for the surface takes the form:

$${{\mathsf{\rho }}_{\mathrm{o}}{\mathrm{K}}_{\mathrm{M}}\frac{\mathsf{\partial }(\mathrm{u},\mathrm{v})}{\mathsf{\partial }\mathrm{z}}|}_{\mathrm{z}=\mathsf{\eta }}=({\mathsf{\tau }}_{\mathrm{wx}},{\mathsf{\tau }}_{\mathrm{wy}})$$

(7)

In Equation (7), η is the free surface elevation in (m) and τ_wx, τ_wy are the wind stress components in Nm⁻² calculated from the equation

τ = C_D ρ_a |W|* W

(8)

where ρ_a = 1.2 kg·m⁻³ is the density of air, C_D is the drag coefficient in Equation (8) computed according to [60]. While (W) is the wind speed component and (K_M) is the vertical mixing coefficient for momentum expressed in m²·s⁻¹.

Daily averaged freshwater discharges were specified for a total of 29 rivers across Ireland, west Britain, and west France. The climatological mean river discharge values were obtained from various sources (

Table 1. Mean annual climatological river discharge values (m³/s) included in the NEA_ROMS.

Country	River Name	Mean Annual Discharge (m3/s)	River Salinity
Ireland	Shannon_Maigue	195.57	0
	Erne	96.72	0
	Corrib	101.87	0
	Barrow-Nore-Suir_Anner	130.87	0
	Boyne	37.63	0
	Blackwater	55.72	0
	Moy	56.16	0
	Slaney	31.83	0
	Liffey	10.93	0
	Lee	27.03	0
	Feale_Galey	26.16	0
	Bandon	15.00	0
Great Britain	Severn_Wye	179.53	0
	Eden_Annan	81.34	0
	Ribble	32.96	0
	Tywi	39.33	0
	Dee	31.03	0
	Lune	35.36	0
	Tamar	22.39	0
	Teifi	28.60	0
	Taw_Torridge	33.74	0
	Derwent	26.01	0
	Leven	43.82	0
	Exe	16.02	0
France	Blavet	29.05	0
	Loire	873.49	0
	Seine	579.11	0
	Dordogne_Garonne	730.72	0
	Adour	44.82	0

), and the location of each river is shown in Figure 1. The NEA_ROMS rivers salinity was set to zero and the river temperatures were not prescribed. Discharge rates for the Irish rivers were obtained from the Office of Public Works (OPW) http://www.opw.ie/ and Electricity Supply Board (ESB) databases. The NEA_ROMS operational river discharge system is evolving. In 2019 our team managed to implement near real time flow rates for Corrib River obtained from OPW. The Corrib River is located inside Galway Bay on the west coast of Ireland. UK rivers data were downloaded from the Centre for Ecology and Hydrology website (https://www.ceh.ac.uk/our-science/projects/national-river-flow-archive). Data from the French rivers were provided from the national database “HYDRO”, run by the French ministry of ecology and solidarity transition. The corresponding model salt flux boundary condition is given according to [61]:

$${{\mathrm{K}}_{\mathrm{H}}\frac{\mathsf{\partial }\mathrm{S}}{\mathsf{\partial }\mathrm{z}}|}_{\mathrm{z}=\mathsf{\eta }}=\mathrm{SSS}(\mathrm{E}-\mathrm{P}-(\frac{\mathrm{R}}{\mathrm{A}}))$$

(9)

where E is the evaporation, P the precipitation, R is the Rivers climatological daily discharge, A is mouth cross-section in m², and SSS in Equation (8) represents the sea surface model salinity.

2.3. Operational Setup

Figure 2 describes the schematic of the NEA_ROMS operational functioning during the simulation period. At the beginning of each year, a new NEA-ROMS run is initialized using interpolated fields from the CMEMS Ocean global model. Tide and river forcing files are created to provide forcing data for the full year ahead. The NEA_ROMS model is then run in hindcast/forecast mode for the year. This entails accessing the CMEMS Ocean weekly analysis run to provide boundary data for the NEA-ROMS hindcast run, while the latest CMEMS Ocean forecast data is used for our forecast boundary forcing. Similarly, the latest ECMWF forecasts are used for atmospheric forcing. The 7-day hindcast is initialized using a restart file from the end of the previous 7-day hindcast. The final output file of the hindcast is used to initialize the first forecast of the current week. A 3-day forecast is created every day using a restart file from the previous forecast to initialize the run. NEA_ROMS state variables (i.e., temperature, salinity, sea surface height, barotropic and baroclinic velocity fields) are saved to netCDF files in hourly snapshots. Selected points in the domain have data saved at the higher frequency of 10 min for use in model validation and also to provide boundary data for two high-resolution coastal models which are also part of the operational system. Model results are published on the Marine Institute website http://www.marine.ie/Home/site-area/data-services/marine-forecasts/marine-forecasts. A rolling month of model output is published to the Marine Institute THREDDS http://milas.marine.ie/thredds/catalog.html and ERDDAP https://erddap.marine.ie/erddap/index.html servers. Monthly mean fields of some parameters (e.g., surface and bottom temperature) are also published to these servers.

2.4. Observations Data and Methods of Analysis

Data used to validate the NEA-ROMS include Sea Surface Temperature (SST) from satellites, tide gauges, ARGO floats, and CTD salinity and temperature profiles.

The satellite SST product used is the European North West Shelf/Iberia Biscay Irish Seas High Resolution ODYSSEA L4 Sea Surface Temperature Analysis dataset which is available for download from the CMEMS website (data product, SST_ATL_SST_L4_NRT_OBSERVATIONS_010_025). The product consists of daily averaged SST values. The data valid for a particular day (from previous day 12:00 to current day 12:00) at a horizontal resolution of 0.02° × 0.02° for a domain covering the European North West Shelf and Iberia-Biscay-Irish Seas. It is a multi-sensor Level 4 analysis which aims to provide an estimate of the night time SST based on original SST observations http://resources.marine.copernicus.eu/documents/QUID/CMEMS-SST-QUID-010-025.pdf. The data from 1 June 2019 were replaced by ODYSSEA product, SST_NWS_SST_L4_NRT_OBSERVATIONS_010_003, https://cmems-resources.cls.fr/documents/PUM/CMEMS-OSI-PUM-010-003-010.pdf.

The authors validated NEA_ROMS Sea Surface Height (SSH) using observed time series from tide gauges around the Irish coast (Irish National Tide Gauge Network) https://data.gov.ie/dataset/irish-national-tide-gauge-network. The sampling frequency of the measured tide data ranged from 5 to 15 min while the model output is produced at a frequency of 10 min. Harmonic analysis was carried out on the observed and modelled time series using the T-TIDE software in MATLAB [62]. In addition, the surge component was calculated following [62], for observed and modelled data at the locations of Irish tide gauges, and the difference statistics are presented. According to [63], the surge (residual) component of sea level is defined as the total water level minus the tide.

Sea level elevation = predicted tide level + storm surge height

(10)

In this study, ARGO float profiles were used to validate the NEA_ROMS temperature and salinity fields. Drifting ARGO floats record vertical profiles of temperature and salinity in the upper 2000 m of the ocean approximately every 10 days and transmit the data via the ARGO satellite system to a number of data processing centers [64,65,66,67]. The data is then made freely available via public access site http://www.argo.ucsd.edu/. Irish Marine Institute (IMI) have used only ARGO data in delayed mode. For the corresponding modelled profile, the authors used the hourly snapshots model output file corresponding to the date the ARGO profile was acquired, and the grid point closest in location to the ARGO profile. The ARGO temperature and salinity were then interpolated onto the modelled vertical profile for that grid point. The authors chose the year 2019 as one that has the best ARGO tracks coverage for the south of the Porcupine Bank, the southern and northern entrance to the Rockall Trough area.

The IMI maintains an archive of all Conductivity Temperature Depth (CTD) measurements taken from the research vessels operated by the Institute. During the period of January 2016 to December 2019, CTDs were acquired from the research vessels in the shelf seas around Ireland with the highest frequency in the months of spring and summer. The authors chose the year 2017 as one that provides the best Irish shelf waters coverage. These data are analyzed to assess the model performance in the shallow, coastal waters around Ireland.

In the study, the authors used the Geostrophic and Ekman Current Observatory GEKCO surface data to validate model velocity fields and Eddy Kinetic Energy (EKE), then compare it with CMEMS Ocean global model (EKE). The data are daily averages 0.25° spatial resolution obtained from Mean Absolute Dynamic Topography (MADT) at the gridded DUACS products. The DUACS system has been producing, as part of the CNES/SALP project, and the Copernicus Marine Environment and Monitoring Service (CMEMS), high quality multi-mission altimetry Sea Level products for oceanographic applications, climate forecasting centers, geophysics and biology communities http://www.aviso.altimetry.fr [68,69]. The geostrophic velocity anomalies presented in this study were deduced from MADT maps with geostrophic approximation as follows:

For NEA_ROMS domain, the mean velocity ${\overline{\mathrm{U}}}_{\mathrm{g}},{\overline{\mathrm{V}}}_{\mathrm{g}}$ in the x and y directions and the anomalies ${\mathrm{U}}_{\mathrm{g}}^{\prime }$, ${\mathrm{V}}_{\mathrm{g}}^{\prime }$ from the mean were computed

$${\mathrm{U}}_{\mathrm{g}}^{\prime } = -\frac{\mathrm{g}}{\mathrm{f}}\frac{\mathsf{\partial }{\mathrm{h}}^{\prime }}{\mathsf{\partial }\mathrm{y}},{\mathrm{V}}_{\mathrm{g}}^{\prime }=\frac{\mathrm{g}}{\mathrm{f}}\frac{\mathsf{\partial }{\mathrm{h}}^{\prime }}{\mathsf{\partial }\mathrm{x}}$$

(11)

where ${\mathrm{h}}^{\prime }$ is the absolute dynamic topography anomaly in meters, g is the gravity acceleration in m/s² and f = 2Ω sin ϕ, where Ω = 2π/T is the earth angular velocity in s⁻¹, T is the earth periodic time = 86,400 s (1 day in seconds), and ϕ is the latitude in degrees. The variance of these velocity anomalies was considered as Eddy Kinetic Energy (EKE), representative of mesoscale activity.

$$\mathrm{EKE} = \frac{ 1}{2} ({\mathrm{U}}_{\mathrm{g}}^{\prime 2}+\prime {\mathrm{V}}_{\mathrm{g}}^{\prime 2})$$

(12)

Model velocities were taken daily at the first sigma layer to compute EKE_{NEA_ROMS} to be exactly like EKE_GEKCO and EKE_CMEMS between 1 January 2017 and 31 December 2017. Velocity time series were subsampled at a daily period of GEKCO MADT as described in [70]. Space and time dependent EKE_GEKCO and EKE_{NEA_ROMS} fields were eventually averaged temporally to produce mean EKE fields.

3. Validation of NEA_ROMS Results against Observations

In this part, the authors will validate and discuss the NEA_ROMS results against available observations described previously in Section 2.4.

3.1. Comparison between NEA_ROMS and ODYSSEA SST during the Simulation Period

Figure 3a–c shows the temporal evolution of the daily mean SST basin averaged statistics for the NEA_ROMS Model-ODYSSEA satellite data for the whole model simulation period from 1 January 2016 until 31 December 2019 in °C. The daily SST bias between the NEA_ROMS and ODYSSEA satellite is shown in Figure 3a. The ODYSSEA SST product is based on night-time data, so for better validation, we rejected all NEA_ROMS daytime SST data. The highest daily bias mean was about 1.8 °C and recorded in mid-June 2018. While the lowest daily bias mean was −0.8 °C in mid-November 2019. The SST for the NEA_ROMS was greater than satellite during the summer months of June and July over the whole simulation period. This may be due to model error as a result of air–sea physics parametrizations during these months. This error has been discussed by [71,72,73]. The lowest bias differences between the model and satellite SST were found from mid-February till the end of April and during October over the whole simulation period (see Figure 3a). The root mean square error (RMSE) between the model and the satellite is described in Figure 3b. This figure suggests that the RMSE between our model SST and the satellite is around 0.5 °C for most of the simulation period. The maximum RMSE was found to be more than 1.5 °C during the summer season for all years except 2019. This may be attributed to problems mentioned above of the model air–sea physics parametrizations during the summer season. Regarding Figure 3c, the authors found that the NEA_ROMS SST is significantly correlated to satellite SST in all simulation period at 95% confidence limit. The highest correlation coefficient was more than 0.97 and recorded in March and April 2018. The minimum correlation coefficient was 0.85 and recorded in mid-June 2018 and December 2019. The annual bias maps between the model and satellite SST are made for the whole simulation period as shown in Figure 4a–d. The annual bias SST maps pattern shows a consistent improvement from the year 2016 until 2019. The maps suggest that minimum bias patterns found in 2019, while the maximum one recorded in 2016. During the simulation period, an effort was made to improve the temperature field by nudging of the surface temperatures towards a reference solution and nudging of the open ocean boundary temperatures and salinity in the relaxation zone towards a reference solution. Referring to Equations (4) and (5), we tuned best values for Γ the relaxation time and heat flux correction terms, ∂Q/∂T to improve our temperature and salinity fields. The NEA_ROMS model appears to overestimate SST by about 1 °C more than ODYSSEA off the European Northwest Shelf as shown in Figure 4a–d but this overestimation pattern has improved in 2019 (see Figure 4d). The authors observed a persistent positive SST bias of more than 1 °C in the Bay of Biscay and along the Iberian Coast except in 2019 as present in Figure 4d. The SST deficiencies in the Bay of Biscay could be due to the model spatial resolution or inaccuracies in the bathymetry data. This region features very deep waters and steep topography. Additionally, the north Atlantic is affected by large surface waves throughout the year, especially during winter time as mentioned by [74]. The surface wave is not included in our model and nonbreaking wave-induced mixing effect could be important in simulating SST as described in [75]. This may lead to significant cooling of the simulated SST [75]. In summary, the spatial distribution of the bias value clearly shows that there are improvements at the model boundary as well as over the whole domain which demonstrates that both changes in the model configuration are contributing to the improvement in model skill.

3.2. Validation of NEA_ROMS Sea Surface Height against Tide Gauge Stations

Observed time series from tide gauges around the Irish coast are compared with model sea surface heights output at the same locations as shown in Figure 5,

Table 2. The amplitudes in meters and phases in degrees for six of the principal tidal constituents calculated, for the measured and modelled data.

Tidal Constituent		Gauge		NEA_ROMS
		Amp	Phase	Amp	Phase
Aranmore	M2	1.22	161	1.29	156
	S2	0.44	194	0.46	189
	N2	0.25	139	0.26	135
	K1	0.12	154	0.14	142
	O1	0.08	359	0.08	5
	Q1	0.03	299	0.03	292
Ballycotton	M2	1.42	148	1.26	146
	S2	0.44	193	0.41	190
	N2	0.26	128	0.23	128
	K1	0.02	178	0.02	178
	O1	0.03	35	0.02	55
	Q1	0.01	340	0.01	348
Ballyglass	M2	1.16	158	1.21	152
	S2	0.41	191	0.43	185
	N2	0.23	136	0.25	131
	K1	0.14	122	0.16	116
	O1	0.09	337	0.08	342
	Q1	0.03	277	0.03	268
Castletownbere	M2	1.12	131	1.09	124
	S2	0.36	161	0.36	155
	N2	0.23	108	0.22	101
	K1	0.04	50	0.05	49
	O1	0.01	279	0.02	277
	Q1	<0.01	185	<0.01	184
Dunmore East	M2	1.38	150	1.20	150
	S2	0.45	199	0.43	198
	N2	0.25	133	0.22	136
	K1	0.04	178	0.04	177
	O1	0.04	29	0.03	49
	Q1	0.01	344	0.01	341
Galway	M2	1.57	141	1.64	137
	S2	0.55	172	0.58	170
	N2	0.32	119	0.34	117
	K1	0.09	76	0.10	79
	O1	0.06	311	0.05	311
	Q1	0.02	261	0.02	243
Howth	M2	1.44	325	1.39	318
	S2	0.41	357	0.41	349
	N2	0.28	297	0.28	289
	K1	0.10	198	0.11	187
	O1	0.08	37	0.07	54
	Q1	0.03	343	0.03	340
Malin Head	M2	1.12	178	1.19	172
	S2	0.42	206	0.44	201
	N2	0.23	156	0.26	151
	K1	0.09	167	0.11	154
	O1	0.07	7	0.07	18
	Q1	0.03	310	0.03	306

and

Table 3. Bias and RMSE of error for modelled minus observed at Irish tide gauge locations in meters.

Tide Gauge Station	NEA_ROMS Bias (m)	NEA_ROMS RMSE (m)
Aranmore	0	0.14
Ballycotton	0	0.16
Ballyglass	−0.01	0.14
Castletownbere	0	0.14
Dunmore East	0.01	0.17
Galway	0	0.15
Howth	0	0.15
Malin Head	0.01	0.15

. The tide gauge stations are distributed around the whole of the Republic of Ireland. Harmonic analysis was carried out on the eight locations using the T-TIDE software in MATLAB [62]. The aim of the analysis was to compare the modelled and measured values of the main tidal constituents (magnitude and phase angle). This analysis demonstrated that the tidal signal in the Sea Surface Height (SSH) data was dominated by three semi-diurnal constituents (M₂, S₂, N₂) and three diurnal constituents (K₁, O₁, and Q₁). Regarding our tidal analysis the authors found that M₂ and S₂ are responsible for most of the tide in our domain in agreement with [76,77,78,79]. The spatial variability in the mean error between model and observations for the two main constituents, M₂ and S₂, is shown in Figure 5a–d. The model shows good skills on the west coast of Ireland and also in the western Irish Sea. The maximum M₂ magnitude error for the west Irish coast was about 0.07 m for Malin Head, Aranmore, and Galway tide gauge stations. However, the model has a higher negative error of about 5% of the total M₂ amplitude in the Celtic Sea recorded in Dunmore East tide gauge station where the M₂ amplitude is consistently underestimated. The S₂ amplitude error shows a broadly similar pattern of M₂. The maximum S₂ error of about (-) 0.03 m (about 6% of the total S₂ amplitude) was recorded in the Celtic Sea at the Ballycotton tide gauge station. One explanation for the poorer skill in the Celtic Sea is an excessive dissipation of tidal energy over the broad expanse of shelf south of Ireland, leading to a reduction in the amplitude of the tidal wave in agreement with [80]. The phase angle errors for both components M₂ and S₂ are negative at most tide gauges stations except in Dunmore East station as shown in Figure 5c,d. The highest negative phase angle error for both components was found to be (-) 7° in the Celtic Sea and also recorded at the Howth gauge station. The aggregated results for the most significant tidal constituents are presented in Table 2. The magnitude difference was very small for all tidal constituents. Ratios of the modelled to the observed amplitudes of the constituents show that the model has good skill at N₂, O₁, and Q₁ though these constituents are much less significant to the total tide than the other three.

Accurate prediction of the surge is vital if a model is to be effective as part of a coastal flooding early warning system. The surge component was calculated following [63], for observed and modelled data at the locations of Irish tide gauges and the different statistics are presented in Table 3. Surge is caused by atmospheric pressure and wind. In order to obtain the surge signal we removed tides. The tidal signal in the Sea Surface Height (SSH) data was dominated by three semi-diurnal constituents (M₂, S₂, N₂,) and three diurnal constituents (K₁, O₁, and Q₁). The bias results were approximately zero for most gauge stations except in Ballyglass, Dunmore East, and Malin Head that were −0.01, 0.01, and 0.01 m respectively. The results show that the RMSE is around 0.15 m for all gauge stations. There are a number of possible contributors to this error, including the slight phase difference in the tidal signal and inaccuracies in the modelled weather fields both in space and time.

3.3. Validation of NEA_ROMS with ARGO Floats Temperature and Salinity Profiles

The locations of ARGO float profiles in the NEA_ROMS domain are shown in Figure 6a–e. The authors used 304 quality checked ARGO profiles data for the year 2019. The bias between modelled and measured ARGO float temperature and salinity profiles are presented in Figure 6a,b. The bias results reveal an agreement between NEA_ROMS and ARGO temperature and salinity profiles for most of the locations. The maximum bias between the NEA_ROMS and Argo floats temperature was 0.8 °C and located at the western open boundary (see Figure 6a). The same locations recorded maximum salinity bias of around 0.15. This may be due to the fact that some of the Argo floats are located too close to the boundary where they are influenced by boundary nudging effects. Additionally, the absence of the wave effect in our model could be an important reason for poor representation of the salinity and temperature across the water column as described in [74]. Biases were used to calculate the RMSE as shown in Figure 6c,d. The maximum RMSE between the NEA_ROMS and Argo floats temperature was 1 °C at the same locations described above. The RMSE magnitude of the temperatures is warmer in the model than measured. While the maximum RMSE for the salinity was about 0.15 and recorded again at the northwestern open boundary edge as seen in Figure 6d. This can be attributed to the reason mentioned above. The NEA model correctly reproduces salinity/temperature profiles in the inner domain and high latitude areas. Figure 6e shows the T-S diagram for both model and ARGO float profiles combined with σ_θ contours, while, the water masses description for the Rockall Trough area according to [24,81] is shown in Figure 7. Figure 6e suggests that ARGO and model output broadly match. The model could describe the high saline dense Mediterranean Water (MW). This water mass is located to the south of the Porcupine Bank at 1000 m depth with a maximum salinity of about 35.8, temperature range from 8 to 10 °C, and a σ_θ of approximately 27.6 in agreement with [82]. In addition, the deep water mass of Labrador Sea Water (LSW) is well developed by the model and ARGO floats T-S diagram. This cold deep water mass is generated over the Labrador basin and slowly makes its way east across the Atlantic Ocean [83]. This water mass is found at 2000 m depth with a salinity of about 34.8, temperature 3 °C, and σ_θ of approximately 27.8 in agreement with [83,84,85] and shown by Figure 6e and Figure 7. The water mass Subarctic Intermediate Water (SAIW) signal is presented by Model and ARGO floats T-S diagram. This water mass is located at the northern entrance to the Rockall Trough at 900 m depth with a maximum salinity of about 35.2 and σ_θ of approximately 27 as shown in Figure 6e and Figure 7, and described by [86]. There was a problem in presenting the surface water mass Eastern North Atlantic Water ENAW at the model T-S diagram. The model does not properly recreate this signal maybe because it mixed with the underneath SAIW water mass. This may be as indirect effect of model excess vertical mixing due to the use of associated parameters with KPP vertical turbulent closure scheme as described in [46,47]. This parameters were used to tune the vertical profiles of currents, temperature, and salinity [46,47]. In conclusion, the comparison with ARGO floats profiles has shown a good capability of the model to reproduce the main structures of the field in the area.

3.4. Validation of NEA_ROMS with CTDs Temperature and Salinity Profiles

To measure the model performance in shelf waters IMI used Conductivity Temperature Depth (CTD) profiles. The method used to compare the model and CTD data is the same as that used for the ARGO data. The authors analyzed 223 CTD profiles From January 2017 to December 2017, CTDs were acquired from the IMI research vessels in the shelf seas around Ireland as shown in Figure 8a–d. This data was analyzed to assess the model performance in the shallow, coastal waters around Ireland. Figure 8a,b presents the bias between modelled and measured CTD temperature and salinity profiles. The temperature bias results propose a good agreement between NEA_ROMS and CTD profiles as shown in Figure 8a. The temperature bias results are almost zero for most locations except one profile situated to the west of the Irish coast. The temperature bias in this single location was found to be 0.9 °C. The salinity bias results suggest an overestimation by the NEA_ROMS of about 0.15 for most profiles located at the southern entrance to the Irish Sea (see Figure 8b). The reason behind this could be related to the lower accuracy of the salinity field for the CMEMS model in this area. The CMEMS salinity field has improved a lot since 2017 as described in [87]. This improvement would obviously have a big impact on our model salinity field. The maximum RMSE between the model and the CTD temperature profiles was found to be 0.8 °C while it was 0.17 for salinity as described in Figure 8c,d. The distribution patterns of the RMSE values clearly show the model robustness in presenting the main water masses at the Irish shelf waters.

The T-S diagram for NEA_ROMS and 223 CTDs profiles describes the Irish coastal waters around Ireland as shown in Figure 8e. The comparison between the CTDs T-S diagram and the model T-S shows a very good agreement. We observed a significant resemblance in the T-S diagram patterns between the CTD data and the model in the whole water column. The NEA_ROMS T-S diagram represents deep water masses well, while it overestimates salinities in some surface shallow coastal waters with maximum centered salinity of about 35.5 and σ_θ of approximately 26.5 as shown in Figure 8b. The reason for the overestimation model salinity is probably a significant underestimation of the discharge of fresh water into the coastal waters of Ireland and west Britain. Another possibility, as previously mentioned, may be due to the use of KPP vertical turbulent closure scheme inside the ROMS model which can affect the vertical mixing as described in [46,47]. This issue is a potential source of deficiencies between model output and observations.

In conclusion, the NEA model correctly reproduces the shelf water masses around Ireland.

3.5. Validation of NEA_ROMS with GEKCO Surface Data

Mean Kinetic Energy (MKE) and EKE calculated from geostrophic current from altimetry GEKCO data and surface current from NEA and CMEMS global models are shown in Figure 9 which was previously described in Section 2.4. The authors combined surface current with MKE maps in Figure 9a–c to validate NEA_ROMS current pattern. The GEKCO data and both models are characterized by the North Atlantic Current and European Slope Current which transports heat and salt from the north-east Atlantic, interacts with the continental shelf slope and forms branches that flow into the North Sea in agreement with [23,35]. Both NEA_ROMS and CMEMS Global models produce similar mean circulation to the GEKCO currents. NEA_ROMS has well-resolved meanders and eddies, which are not well-resolved by CMEMS global (see Figure 9b). For example, the eddies and meanders associated with the North Atlantic Current at the middle of NEA_ROMS domain seem to be realistic as shown before by [32,35,88]. The NEA_ROMS has a higher resolution and more detailed bathymetry in this area than the CMEMS model, and this is probably the reason that it is better able to resolve the eddy field in the area. Figure 9a–c suggests that CMEMS global model generally has the largest MKE field, except for the North Ireland region. However, CMEMS global assimilates observations such as altimetric sea level, temperature and salinity profiles, and corrects toward satellite SST at the surface, as described in [87], so that MKE could be larger at times due to the correction of physics uncertainties [73,89]. Figure 9e–f represents the mean EKE fields for GEKCO data, NEA_ROMS, and CMEMS global, respectively. The basin averaged EKE for GEKCO data, NEA_ROMS, and CMEMS were 0.0063, 0.0153, and 0.0147 m² s⁻², respectively. The EKE comparison suggests the NEA_ROMS has the highest EKE field. The finer model resolution and the better river discharges of NEA_ROMS may produce this change according to [73,90]. Overall, NEA_ROMS reproduces the main regional circulation patterns, improving the dynamics forced by finer resolution and more detailed bathymetry in this area.

4. The Irish Coast Current (ICC)

Whilst it is fundamental that a hydrodynamic model is quantitatively validated, it is also necessary that it reproduces the important hydrographical features within its domain. To that end, an assessment is made of the model’s ability to reproduce the summer Irish Coastal Current (ICC) [10,11,91]. The ICC flows northward parallel to the western coast of Ireland [92]. The flow of ICC depends on the local wind forcing, so the Atlantic Water then either enters the Irish Sea through the North Channel or, more likely, continues flowing northward producing the ICC until the Hebrides [92]. Additionally, the ICC shows a shelf-located meandering surface current of 10–20 cm·s⁻¹ which, has a full depth presence over the 100–120 m contours [93,94]. The authors followed the methodology of [11] to produce the annual maps of the summer ICC simply by averaging two months (July and August) of the mid-depth residual NEA_ROMS velocity field for the Celtic Sea and Irish Shelf. We used the model output hourly snapshots velocity field. To produce annual maps, we get velocities away from ICC region out of the picture.

Figure 10 displays the annual results for this analysis over the whole simulation period from 2016 until 2019. The velocity of the current is found to be in the range of 6–20 cm s⁻¹ which is in agreement with reported values in the published literature [10,11]. The localized “hotspots” of increased velocities near headlands on the western coast are presumably induced by the tidal current while the weakening of the current north of the Shannon may be indicative of tidal mixing being more dominant than stratification [11,93]. The map patterns suggest a well-developed ICC for all years except for 2016. The highest velocity for the ICC was about 19 cm s⁻¹ and recorded in 2017, 2018, and 2019 as shown in Figure 10b–d. There is an interesting difference in 2019 where the ICC is very strong and well defined along most of the western Irish coast (see Figure 10d). In conclusion, the results show that the NEA_ROMS has reproduced the ICC front for the whole simulation period.

5. Conclusions and Further Discussion

This paper presents an operational model for Irish waters NEA_ROMS. The model is based on ROMS and is forced by operationally available atmospheric and boundary data. The system is robust following several years in operation and model data is made freely available via THREDDS server. The authors validated and calibrated the NEA_ROMS with available observations in our domain area. The observations were the SST obtained from satellite covering the simulation period, time series of eight tide gauges around the Irish coast obtained from the Irish National Tide Gauge Network, Argo, and CTD temperature and salinity profiles.

The SST for the NEA_ROMS bias was greater than the satellite during the summer months of June and July over the whole simulation period. The highest daily bias mean was about 1.8 °C and recorded in mid-June 2018. While the lowest daily bias mean was −0.8 °C in mid-November 2019. The authors attribute this to the model error as a result of air–sea physics parametrizations as discussed by [71,72,73]. The highest correlation coefficient was found to be more than 0.97 and recorded in March and April 2018. While, the minimum correlation coefficient was 0.85 and recorded in mid-June 2018 and December 2019.

The model has a higher negative error of about 5% of the total M₂ amplitude in the Celtic Sea recorded in Dunmore East tide gauge station. The S₂ amplitude error showed a broadly similar pattern to M₂. The maximum S₂ error was about 6% of the total S₂ amplitude recorded in the Celtic Sea at the Ballycotton tide gauge station. That was due to the poorer skill in the Celtic Sea due to excessive dissipation of tidal energy over the broad expanse of shelf south of Ireland, leading to a reduction in the amplitude of the tidal wave in agreement with [80].

The validation of NEA_ROMS with ARGO floats revealed the model skills in representing the main water masses in particular MW and LSW for the area located to the south of the Porcupine Bank and the entrance to the Rockall Trough. The location of some Argo floats too close to the boundary has increased the temperature and salinity RMSE because they may have been influenced by model boundary nudging effects.

The maximum RMSE between the model and the CTD temperature profiles was found to be 0.8 °C while it was 0.17 for salinity. The NEA model correctly reproduces the shelf water masses around Ireland.

The use of KPP vertical turbulent closure scheme in the ROMS model can affect the vertical mixing of the salinity and temperature profiles as described in [46,47]. This issue was a potential source of deficiencies between model predictions and observations.

The meanders and eddies at the model domain were well resolved by NEA_ROMS which were not well resolved by CMEMS global. The NEA_ROMS has a higher resolution and more detailed bathymetry in this area, which may be the reasons why the model resolved these meanders and eddies well.

In 2019 the ICC was very strong and well defined along most of the western Irish coast. In conclusion, the results showed that the NEA_ROMS well reproduced the ICC front for the whole simulation period.

Future work requires the use of high frequency near real time river discharges. The inclusion of data assimilation in the NEA_ROMS operational system is considered to be the next logical step in its ongoing development. Recent applications also include the modelling studies of the biogeochemical cycling. Additionally, indices supporting the implementation of the Marine Strategy Framework Directive (MSFD) are usually obtained from in-situ data. This is a major difficulty in oceanic areas where data are scarce. Validated numerical models can fill this gap. Tools will be developed to obtain indices of interest directly from model results, e.g., areas of upwelling, fronts, eddy index, primary production, trophic status, or even conditions for propagation of noise.