**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–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<sup>3</sup> s<sup>−</sup><sup>1</sup> [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].

**Figure 1.** Bathymetric map of the area covered by the model. Major rivers included in the model are marked by red circles.

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 grea<sup>t</sup> 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–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–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–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 di ffusion coe fficients. Vertical di ffusion 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–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 rx0 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{\boldsymbol{\wp}\_{\text{nested}}} = \left[ \sqrt{\frac{\mathbf{g}}{\mathbf{H}\_{\text{nested}}}} (\eta\_{\text{nested}} - \eta\_{\text{nestic}}) \right] - \overline{\boldsymbol{\wp}\_{\text{nested}}} \frac{\mathbf{H}\_{\text{nestic}} + \eta\_{\text{nestic}}}{\mathbf{H}\_{\text{nested}} + \eta\_{\text{nested}}} \tag{1}$$

where ϕn<sup>=</sup> 1 H+η η H ϕndz is the normal barotropic velocity at the open boundaries and g the gravity, whiletheinflowwillbe

 at ϕnnesting = ϕnnested

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

$$\int\_{n\_2}^{n\_1} \int\_{\text{nested}}^{n\_{\text{nested}}} \varphi\_{\text{nned}} \, \text{dzdn} = \int\_{n\_2}^{n\_1} \int\_{\text{nension}}^{n\_{\text{nossing}}} \varphi\_{\text{nensing}} \, \text{dzdn} \tag{2}$$

where ϕ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{\partial \mathbf{T}\_{\text{nested}}}{\partial \mathbf{t}} + \boldsymbol{\varphi}\_{\text{nreset}} \frac{\partial \mathbf{T}\_{\text{nested}}}{\partial \mathbf{y}} = 0; \frac{\partial \mathbf{S}\_{\text{nested}}}{\partial \mathbf{t}} + \boldsymbol{\varphi}\_{\text{nreset}} \frac{\partial \mathbf{S}\_{\text{nested}}}{\partial \mathbf{y}} = 0 \tag{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 ϕ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{\partial \gamma}{\partial t} = \text{r.h.s.} \frac{1}{\Gamma} \left( \gamma - \gamma^{\text{CMEMS}} \right) \tag{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 (M2, S2, N2, K2, K1, O1, P1, Q1, Mf, 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

$$\rho\_{\rm o} \mathbf{C}\_{\rm P} \mathbf{K}\_{\rm H} \frac{\partial \mathbf{T}}{\partial \mathbf{z}}\Big|\_{\mathbf{z}=\eta} = \mathbf{Q}\_{\rm T} + \frac{\partial \mathbf{Q}}{\partial \mathbf{T}} \left( \text{SST}\_{\rm CEMS} - \text{SST}\_{\rm m} \right) \tag{5}$$

where Q**T** is the total net heat flux according to the formula:

$$\mathbf{Q\_T = Q\_S - Q\_b - Q\_e - Q\_h} \tag{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<sup>−</sup>2· ◦C−<sup>1</sup> for NEA-ROMS. SST m is the model sea surface temperature and salinity and SSTCMEMS is the reference SST. Q**<sup>S</sup>**, Q**<sup>e</sup>**, Q**<sup>b</sup>**, and Q**h** are the short, latent, back, and sensible heat fluxes in <sup>W</sup>/m2. (K H) is the vertical mixing coe fficients for tracers in (m<sup>2</sup>·s<sup>−</sup>1). 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:

$$\rho\_{\rm o} \mathbb{K}\_{\rm M} \frac{\partial(\mathbf{u}, \mathbf{v})}{\partial \mathbf{z}} \bigg|\_{\mathbf{z} = \boldsymbol{\eta}} = (\boldsymbol{\pi}\_{\rm W \mathbf{x}, \prime} \boldsymbol{\pi}\_{\rm W \mathbf{y}}) \tag{7}$$

In Equation (7), η is the free surface elevation in (m) and τwx, <sup>τ</sup>wy are the wind stress components in Nm−<sup>2</sup> calculated from the equation

$$
\pi = \mathcal{C}\_D \not\models \mathbb{W} | \mathbb{W} | \prime \; \mathcal{W} \tag{8}
$$

where ρa = 1.2 kg·m<sup>−</sup><sup>3</sup> is the density of air, CD is the drag coefficient in Equation (8) computed according to [60]. While (W) is the wind speed component and (KM) is the vertical mixing coefficient for momentum expressed in <sup>m</sup>2·s<sup>−</sup>1.

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), 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/nationalriver-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]:

$$\left. \mathbf{K}\_{\rm H} \frac{\partial \mathbf{S}}{\partial \mathbf{z}} \right|\_{\mathbf{z} = \eta} = \text{SSS}(\mathbf{E} - \mathbf{P} - \left( \frac{\mathbf{R}}{\mathbf{A}} \right)) \tag{9}$$

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


**Table 1.** Mean annual climatological river discharge values (m<sup>3</sup>/s) included in the NEA\_ROMS.
