High-Resolution Model of Clew Bay – Model Set-up and Validation Results

1. Introduction

Clew Bay is a large bay with an area of about 176 km² (i.e., 16 km from east to west and 11 km from north to south) situated on the west coast of Ireland, characterized by a large number of islands (see Fig. 1). The bay is bounded by Croagh Patrick Mountains to the south and Achill Island to the north (Fig. 1). The islands in Clew Bay are partially drowned drumlins, i.e., elongated, steeply sloping hills that were sometimes called “whalebacks” [1]. They were formed when glaciers reshaped the landscape during the last ice age [1]. Several of the hills on the mainland around the bay are similar drumlins [2]. These glacial formations vary considerably in size, ranging from large islands on which residences and pastures stand to small mounds on the seafloor [1]. The numerous islands give rise to shallow straits and lagoons through which deep channels flow [2]. Erosion of existing and submerged drumlins with their coarse glacial deposits gives rise to a heterogeneous sediment environment [2,3]. Clew Bay has three hundred and sixty-five islands, the largest one, Clare Island, protects the entrance to the sheltered bay [4], as shown in Figure 1.

Clew Bay is an important national region responsible for aquaculture (i.e., finfish farming, oysters farming, mussels, and others) industry in Ireland [5]. In 2020, there were 22 aquaculture-related businesses in Clew Bay [6], the vast majority of which were involved in oyster farming, although shellfish and finfish farming are also represented locally https://bim.ie/wp-content/uploads/2022/05/Clew-Bay-Report-SPREADS.pdf. Finfish and oyster farming extends from the Westport River to the Burrishoole Fishery (Loughs Feagh and Furnace) near Newport and is an important national region [6] (Figure 1).

There is a need for hindcasts and forecasts of hydrodynamic properties for the Clew Bay region by a wide range of stakeholders, e.g., scientists working on salmon migration. To date, there are no published modeling studies for Clew Bay.

As regards other recent modelling studies of Irish waters, the set-up and validation results of a high-resolution numerical operational model developed at the Irish Marine Institute covering the northeast Atlantic with emphasis on Irish waters were presented in [7]. The model is based on ROMS with a horizontal resolution of 1 km and runs operationally. The meanders and eddies in the model domain were well resolved. The authors attribute this to the high resolution of the model and detailed bathymetry in the study area. In addition, [8] has presented a high-resolution (1 km) 3-D ocean model based on ROMS for southwestern Irish waters. The simulation of the Irish hindcast model was validated and calibrated with available observations in the geographic domain of the model. The model has been shown to be robust when compared to CTD vertical temperature and salinity profiles.

Furthermore, in [9], a 3D model of the coastal ocean was created using simple equations for Bantry Bay. The model was developed, validated, and implemented operationally.

The authors have found some other modeling studies of the Irish shelf waters, e.g. [10,11,12,13,14,15,16,17,18,19]. Moreover, in [20] preliminary results of a high-resolution 3D numerical simulation with a wetting and drying scheme based on the Regional Ocean Modelling System (ROMS) were used to study the tidal circulation in Kilmakiloge Harbour. This bay is located on the southern shore of Kenmare Harbour. The model showed good skill and a high correlation coefficient with respect to the available observations, especially for temperature and tides. However, the model overestimated mixing in Kilmakilloge Harbour.

To the best of the authors knowledge, the presented model is the first 3D numerical model developed for Clew Bay. This study reports the validation results of the preoperational model. The model was initialized on June 7, 2017, and ran until January 1, 2019. The model has a horizontal resolution of 80 m and 15 vertical sigma levels. The model was validated using currents recorded by an Acoustic Doppler Current Profiler (ADCP), temperature and salinity profiles from the Irish Environmental Protection Agency (EPA) Ireland website, and water level records from a tide gauge. The Clew Bay net flow and residence time were estimated. Section 2 describes the model implementation and nesting procedures. Section 3 presents the validation against observational data, describes the general circulation in the Clew Bay region, estimates of the net flow and residence time, and discusses the results before presenting the conclusions in the last section.

2. Model Design, Description and Implementation

The model is based on version 3.7 of the code Regional Ocean Modelling System [21]. All model equations are written in rectangular coordinates and include spatially and temporally varying horizontal eddy viscosity and diffusion coefficients [22]. The Clew Bay model has a horizontal grid with a resolution of 80x80 meters and 15 vertical sigma layers. The high-resolution bathymetry of the Clew Bay model is from Ireland's Integrated Mapping for the Sustainable Development of Ireland's Marine Resource (INFOMAR) database (www.infomar.ie; downloaded on June 19, 2022) and the European Marine Observation and Data Network (EMODnet) bathymetric dataset (https://www.emodnet-bathymetry.eu/). Minimal smoothing of the bathymetry was performed using a linear programming method [23]. The model also has a wetting and drying scheme according to [24]. The scheme specifies that some points of the ocean grid may be completely “wet” or “dry”. These “dry” points may not have zero depth; instead, they are covered by a thin film of water so that the equations of motion can be calculated at all grid points [25]. The direction of flow is determined at each of the faces of a tracer cell, and the flow and velocity through the face are set to zero if the depth of the upstream tracer cell is less than or equal to a user-specified critical depth D_crit [24,25]. The authors set the D_crit to 0.25 m. The turbulence mixing scheme of the model is a k-ε parameterization implemented by the Generic Length Scale (GLS) scheme [26,27]. An upstream third-order scheme with implicit mixing was used for the horizontal advection of momentum [28], whilst multidimensional positive definite advection transport algorithm (MPDATA) was used for horizontal and vertical advection of tracers [29]. Bottom stress is applied using the logarithmic “wall law” with a constant roughness length of 0.01 m.

The Marine Institute's operational model for the Northeast Atlantic (NEA_ROMS) provides the initial conditions for the Clew Bay model of the three open boundaries [7]. The initial and boundary parameters of the model include temperature, salinity, baroclinic and barotropic velocity components, and sea surface height. The temporal resolution of the boundary conditions is 10 minutes and includes the tidal signal. The nesting was designed so that the volume transport across the open boundary of the nested model matches the volume transport of the NEA _ROMS model; this technique was described in [7].

The atmospheric data for the computation of the surface forcing are taken from the global high resolution (0.125^o) atmospheric model run by the European Centre for Medium Range Weather Forecasts (ECMWF) at a 3 -hour frequency. The atmospheric fields used are air temperature, relative humidity, wind speed at 10 m above sea level, mean sea level pressure, cloud cover, total precipitation, and solar radiation (surface and net longwave). Wind stress, heat fluxes, and evaporation rates are calculated using an interactive mass formula that uses atmospheric data, as described in [8]. Daily averaged freshwater discharges were specified for eight major freshwater sources (see Table 1). Data for four rivers were obtained from the European Hydrological Forecasts for the Environment (Swedish Meteorological and Hydrological Institute (E-HYPE SMHI)); https://hypeweb.smhi.se/explore-water/historical-data/europe-time-series/) [30], while the remaining rivers were obtained from the Irish Environmental Protection Agency (EPA) Ireland website. . Daily climatological flows for the four rivers in the E- HYPE (Owenwee, Burrishole, Owengrave, and Mayour) were calculated from available flow time series over a 30-year period (1981 - 2010). Daily flow rates for the Newport, Westport, Owennabrockagh, and Bunowen rivers were determined from time series from EPA for a 10-year period (2006-2016). Salinity of incoming freshwater was set to zero, and temperature was calculated from monthly temperatures for the period 2000-2010, available at HYPE-SMHI: (https://hypeweb.smhi.se/explore-water/historical-data/europe-time-series/).

The model was initialized on June 7, 2017 and ran until January 1, 2019. The output consists of temperature, salinity, sea surface height, barotropic, and baroclinic velocity fields and is stored in netCDF files in hourly snapshots. In addition, the output is stored at a frequency of 10 minutes at selected locations in the area for use in model validation.

Data used to validate the 1.5 year hindcast simulation (i.e., June 2017–1^st January 2019), included two ADCPs, four EPA temperature and salinity profiles and Roonagh Tide Gauge Station (see Figure 1).

ADCP is a hydro acoustic current measurement device anchored near the bottom, similar to sonar, and used to measure current velocities over a range of depths, using the Doppler effect of sound waves backscattered from particles in the water column [31]. Two ADCPs were deployed from July 25 to December 20, 2017, and measured currents at a frequency of 12 minutes. The locations of the ADCPs were north of Clare Island at 9.97°W, 53.84°N, where the local water depth is around 40 meters, and in the inner bay at 9.67°W, 53.86°N, where the local water depth is about 15 meters, see Figure (1). The two ADCPs also measured water temperature at the seafloor. Validation was performed for the temperature and barotropic velocity components over the above time period.

Comparison with the barotropic velocity components of the ADCPs rather than the absolute velocity components of the ADCPs reduces the uncertainties in the ADCP velocities, as described in [32,33,34]. We calculated the ADCPs barotropic velocity components ($\overline{\mathrm{U}},\overline{\mathrm{V}}$) by integrating the ADCP velocity components (U, V) related to the water column resolved by the ADCP (H) in (meters) using the following equation:

$$(\overline{\mathrm{U}},\overline{\mathrm{V}}) =\frac{1}{\mathrm{H}}{\displaystyle \underset{-H}{\overset{0}{\int }}(\mathrm{U},\mathrm{V}) \mathrm{d}\mathrm{z}}$$

(1)

Roonagh Tide Gauge Station data for water levels are from 14^th August 2017 to 2^nd January 2018 and were obtained from the Irish National Tide Gauge Network https://data.gov.ie/dataset/irish-national-tide-gauge-network. This tide gauge station is located at Roonagh Pier in the south of the island of Clare and recorded water levels at 6-minute intervals from 14^th August 2017 to 2^nd January 2018. Harmonic analysis of the observed and modelled sea surface heigh time series was performed using T-TIDE software in MATLAB [35]. The objective of the analysis is to compare the modelled and measured values of the main tidal components (magnitude and phase angle).

Data for the temperature and salinity profiles are collected as part of the EPA national water quality monitoring programme (https://www.epa.ie/publications/monitoring--assessment/freshwater--marine/irelands-national-water-framework-directive-monitoring-programme-2019-2021.php). Data were collected using vertical profiles at each station. A data-sonde was lowered from the surface to the bottom of the water column and each parameter is recorded at ~1m intervals where total depth was less than 10m and ~2 m where the total depth was >10m. The data-sondes used were Hydrolab DS5 or Hydrolab HL7 sondes. Each instrument was equipment with sensors to record depth (m), salinity, pH, optical dissolved oxygen (% saturation), turbidity (NTU) and chlorophyll A(µg/l). Data was recorded to a Hydrolab Surveyor handheld datalogger. Salinity measurements were calibrated against KCL standards of known conductivity. The EPA station locations are listed in Table 2.

The Taylor diagram [36] was used to show the validation of the model results. In the Taylor diagram, the correlation coefficient (R), the centered mean square error (RMSD), and the standard deviation (σ_n) are represented by the location of a point representing the test plot relative to the reference point (i.e., Observations). The radial distance of the test point from the origin gives σ_n, its azimuthal position gives R, and its distance from the reference point gives the centered RMSD [36].

3. Results and Discussion

3.1. Validation of the Clew Bay Model against Observations

In this part, we discuss the validation results of the Clew Bay hindcast simulation with available ocean observations described in Section 2.1.

3.1.1. Validation with ADCPs

Figure 2 shows the comparison of the time series between the Clew Bay model and ADCPs bottom temperature in [^oC] for the period from 25^th July to 20^th December 2017. In August and September, the model bottom temperature was about 1.3 ^o C warmer than the ADCP in the inner bay, as shown in Figure 2a. The smallest differences (< 0.2 ^oC) between the Clew Bay model and the inner bay ADCP bottom temperature were observed for the remainder of the Inner Bay ADCP period (i.e., September through December 2017; (Figure 2a)). This period shows clear agreement between model bottom temperature and observation within Clew Bay. Comparison between the Clew Bay model and the Clare Island ADCP temperature showed that the model is in good agreement with the ADCP from September 21^st to 25^th December 2017 (i.e., differences are less than 0.2 °C), as shown in Figure 2b. The largest temperature differences (> 1 ^oC) between the Clew Bay model and the Clare Island ADCP were observed between September 10 and 20, 2017 (see Figure 2b). During this period, the model was cooler than the Clare Island ADCP.

Figure 3 shows the comparison between the Clew Bay model and ADCP bottom temperature in the form of a Taylor diagram and is based on a 5-month period spanning from 25^th July 2017 to 25^th December 2017 with a frequency of 10 minutes. The Taylor diagram generally shows a very high correlation (i.e., R > 0.95) between the model and ADCP bottom temperature. Figure 3.a shows that the correlation coefficient and centred Root Mean Square Difference (RMSD) between the model and ADCP water temperature in the inner bay are nearly 0.99 and 0.5, respectively. For the Clare Island ADCP (Fig. 3b), the correlation coefficient was 0.97 and the RMSD was nearly 0.53. The Taylor diagram also shows that the simulated bottom temperature variations are close to those observed. This is due to the correct model forcing and the fine horizontal (80 m) and vertical resolution (15 layers) of the Clew Bay model, which help to resolve the temperature in the model domain accurately.

Figure 4a–d presents Taylor diagrams comparing the barotropic velocity components (u, v) between the Clew Bay model and the ADCPs in Inner Bay and Clare Island, with the corresponding correlation coefficients and RMSD. The model results show that the model is significantly correlated with the barotropic velocity components (u, v) of the ADCPs over the above time period with a confidence limit of 95%. The correlation coefficients and RMSD between the model and the ADCPs in Inner Bay and Clare Island for the eastern u-components were [0.85, 0.14] and [0.81, 0.12], respectively, as shown in Figure. 4a,b. For the northern v-components, the correlation coefficients and RMSD between the model and ADCPs are [0.25, 0.03] and [0.82, 0.07], respectively, (see, Figure 4b). The barotropic velocity components (u,v) of the model for the Clare Island site are closer to the observation than those for the inner bay (Figure 4a,b). This can be attributed to the effect of coastal waves; this feature is not implemented in our model [37]. This coastal wave effect is smaller in the Clare Island area with a depth of about 40 m, but larger in shallow areas such as the inner bay with a depth of 15 m [38]. In general, the highest correlation coefficients were found for the u component between model and ADCPs. The modelled meridional current velocity (v-component) in the inner bay is much weaker than the observed value, which could be the reason for the (poor) correlation of the v-component values.

Overall, validation of the Clew Bay model with the ADCPs during the deployment period showed good agreement with both temperature and barotropic velocity components.

3.1.2. Validation of the Model with In-Situ Temperature and Salinity Vertical Profiles

In this part, we discuss the validation results of the Clew Bay simulation with the available temperature and salinity profiles obtained from the EPA and recorded on 20 September 2017 and described in Section 2. Figure 5a–d depicts the validation of simulated temperature, salinity, and density for the profiles at four stations CW030, CW110, CW130, and CW140. The simulated temperature, salinity, and density profiles are similar to the observations (see Figure 5). The RMSE between model and observation for temperature varies between 0.24 and 0.5 ^oC. The maximum temperature RMSE is 0.5 ^oC at station CW110 (Figure 1 and Figure 5b), located in the north of the inner bay, while the minimum RMSE is at a shallow (i.e., depth ~ 9 m) station CW030 in the south of the inner bay (Figure 1 and Figure 5a). The vertical structure of salinity in Clew Bay is generally well reproduced by the model (Figure 6). The salinity at the shallow station CW030 is almost homogeneous and shows the same values. The minimum RMSE of salinity (~0.08) is observed at station CW130 almost in the middle of the inner bay (Figure 1 and Figure 5c). The similarity of the structure between the model and the observed salinity indicates the correct implementation of freshwater in the model. Model density agrees with observations at all stations, especially in the upper 10 meters, and RMSE varies from 0.06 to 0.23. In summary, the analysis shows that the Clew Bay model is well able to reproduce vertical profiles in the study area.

3.1.3. Validation with Roonagh Tide Gauge Station

Figure 6 demonstrates the statistical comparison between the model and the Roonagh tide gauge station based on the Taylor diagram for the period from August 14, 2017, to January 2, 2018. We found that the model is significantly correlated (95% confidence level) with the SSH value of the tide gauge over the above period across the Taylor diagram. The correlation coefficient between the model and the tidal SSH level was 0.99 with an RMSD of 0.09 m. The Taylor diagram also shows that the amplitude of the simulated SSH fluctuations is close to that observed. These results highlight the robustness of the SSH model for Clew Bay, which is probably due to the better boundary condition for the water flow in the model mentioned in [39].

This analysis showed that the tidal signal in the SSH data is dominated by three semidiurnal components with three diurnal components. The components are M₂, S₂, N₂, K₁ O₁, and Q₁. Table 3 shows the comparison between the amplitude and phase angle of the modelled and measured main tidal components with the differences between them (model observation). This analysis showed that M₂ and S₂ are responsible for most of the tides in the area, in agreement with [7,40,41,42,43]. The magnitude differences between the Clew Bay model and the Roonagh tide gauge were very small for all tidal components. The amplitude difference varied from 0.05 m for M₂ (i.e., less than 4% of the total M₂) to 0 for Q₁. There were some small differences in phase angle for the K₁ and O₁ diurnal components. The highest phase angle difference was [+ 9.19◦] for component K₁, while the lowest was [- 0.43◦] for M₂ (see Table 3). The differences between the modelled and observed amplitudes of the components show that the model has good agreement for M₂, S₂, N₂, and Q₁. In conclusion, a very good agreement was obtained for all tidal components.

3.2. Clew Bay Current Patterns

Here, we describe the monthly averaged barotropic current patterns in the Clew Bay model during winter, represented by January, and summer, represented by July. In addition, we present the residual barotropic current on September 10, 2018 (spring tide) and October 2, 2018 (neap tide). The residual current is the current with the tidal signal removed and it was approximated here by averaging 1 hourly current over 25 hours [44,45]. The residual current is the result of several processes such as density-driven current and wind-driven current with tidal signal largely removed.

Figure 7a,b shows the monthly water barotropic velocity fields for the Clew Bay model in winter (January) and summer (July). There are similarities found between all months in the flow patterns of the Clew Bay main circulation. Water flows into Clew Bay south of Clare Island, while north of Clare Island the flow is outward with a cyclonic circulation inside the bay. The current entering the bay is deflected to the right due to Coriolis force. The flow south and north of Clare Island toward and out of Clew Bay is relatively stronger in winter (> 0.1 m/sec) (Figure 7b). This current south and north of Clare Island is weaker in summer (< 0.1 m/sec) because the wind has less influence. The northward flow north of Clare Island is bounded in both seasons by the extent of a small cyclonic (counterclockwise) circulation region occupying geographic positions 53° 51' N and 10° 00' W.

Figure 8a,b depicts residual currents at spring tide on September 10, 2018 and at neap tide on October 2, 2018. Residual barotropic currents at neap tide in Clew Bay are relatively small (< 0.1 m/sec) compared to residual barotropic spring tide (> 0.1 m/sec). The maximum residual barotropic current velocity of about 0.2 [m/sec] is located north of Clare Island at a latitude of ~ 53° 50' N and a longitude of ~ 10° 00' W (Figure 8a). During the neap tide, a slightly stronger longshore current (~0.1 m/sec) flows along the southern coast of Clew Bay (Figure 8b). Overall, the good representation of the barotropic ocean currents is due to the high model resolution (80 m), which helped to identify the spatial variation and resolve the small-scale gradient well.

3.3. Estimation of the Net Flow through Clew Bay

One of the objectives of our study is to estimate the rate of inflow across Clew Bay and the residence time for different months using our model results from the previous section. This information is of great importance to a variety of stakeholders, especially scientists working on salmon migration. Figure 9a,b shows the average monthly water inflow through the Clew Bay cross-section across the southern channel. The inflow rate is estimated from the barotropic velocity of the model times the cross-section area (i.e., depth x dy) over the line shown in Figure 9 map. The maximum water inflow (> 6x10³ m³/s) is observed in January, while the minimum (< 1 x10³ m³/s) occurs in March 2018. There is an interesting coincidence of the maximum and minimum values of wind speed with the maximum and minimum values of inflow rate in the cross section of the bay (Figure 9a,b). The prevailing direction for all months except June and March is southwest (i.e., 270◦ > direction > 180◦; Figure 9b). The south-westerly wind increases the inflow of water into the bay. While it is northwesterly in March (i.e. 270◦ > direction > 360◦) and southeasterly in June (i.e. 180 >direction > 90). The water inflow rate in the bay is strongly dependent on the wind speed and direction as demonstrated by Figure 9 a,b.

4. Conclusions

In this paper, we present for the first time a preliminary result of a high-resolution one-way nested hindcast simulation based on ROMS for the Clew Bay region. The model is driven by lateral boundary conditions taken every 10-minutes from the NEA _ROMS model [7] and atmospheric forcing 3-hourly ECMWF surface fields. Eight freshwater sources were specified, and a wetting and drying scheme implemented. The simulation of the Clew Bay model was validated and calibrated with available real observations (e.g., ADCP, vertical salinity and temperature profiles, and tide gauges) in the geographic area of the model domain.

The correlation coefficient and (RMSD) between the model and ADCP bottom water temperature in the inner bay were 0.99 and 0.5 ^oC, respectively. For the ADCP at Clare Island bottom temperature, the correlation coefficient was 0.97 and the RMSD was nearly 0.53 ^oC. The Taylor diagram also showed that the simulated temperature variations were close to those observed. We attribute this to the correct model forcing and the fine horizontal resolution (80 m). The Taylor diagrams for the barotropic velocity components (u,v) showed that the model for the Clare Island ADCP site is closer to the observations than the model for the Inner Bay ADCP. This could be due to the effect of coastal waves; this feature is not implemented in Clew Bay model [37]. This coastal wave effect according to [38] is smaller in the deep area (i.e., Clare Island area with a depth of about 40 m), but larger in shallow areas such as the inner bay with a depth of 15 m. In addition, the highest correlation coefficients (R >0.80) for the u component between model and ADCPs were found in Clare Island and Inner Bay. The modelled meridional current velocity (v-component) in the inner bay is much weaker than the observed value, which could be the reason for the (poor) correlation (R~0.2) of the v-component values. The resemblance in structure between the model and the observed EPA salinity profiles indicates that freshwater was correctly implemented in the Clew Bay model. The model density agreed with observations at all stations, particularly in the upper 10 meters, and the RMSE ranged from 0.06 to 0.23. The Clew Bay model was able to reproduce vertical profiles in the study area when compared to the EPA. Moreover, the Taylor diagram revealed that the amplitude of the simulated SSH fluctuations is similar to that observed from Roonagh tide gauge station. The model's correlation coefficient with the tidal SSH level was 0.99 with an RMSD of 0.09 m.

There were similarities found between all months in the circulation patterns of Clew Bay. Water flows into Clew Bay south of Clare Island, while north of Clare Island flows out with a cyclonic motion inside the bay. The current south and north of Clare Island toward and out of Clew Bay was stronger in winter. The northward flow north of Clare Island was bounded in winter and summer 2018, by the extent of a small cyclonic (counterclockwise) circulation region near 53° 51’ N and 10° 00’ W.

The maximum water inflow (> 6x10³ m³/s) was found in January, while the minimum (< 1 x10³ m³/s) occurred in March 2018. We noticed a match between the maximum and minimum values of wind speed with the maximum and minimum values of inflow rate. During the winter of 2018, when a strong southwest wind predominated, the inflow of water into the bay was at its highest. The bay's water inflow rate was highly influenced by the wind's velocity and direction.