Comparative Study on Assimilating Remote Sensing High Frequency Radar Surface Currents at an Atlantic Marine Renewable Energy Test Site

Abstract

A variety of data assimilation approaches have been applied to enhance modelling capability and accuracy using observations from different sources. The algorithms have varying degrees of complexity of implementation, and they improve model results with varying degrees of success. Very little work has been carried out on comparing the implementation of different data assimilation algorithms using High Frequency radar (HFR) data into models of complex inshore waters strongly influenced by both tides and wind dynamics, such as Galway Bay. This research entailed implementing four different data assimilation algorithms: Direct Insertion (DI), Optimal Interpolation (OI), Nudging and indirect data assimilation via correcting model forcing into a three-dimensional hydrodynamic model and carrying out detailed comparisons of model performances. This work will allow researchers to directly compare four of the most common data assimilation algorithms being used in operational coastal hydrodynamics. The suitability of practical data assimilation algorithms for hindcasting and forecasting in shallow coastal waters subjected to alternate wetting and drying using data collected from radars was assessed. Results indicated that a forecasting system of surface currents based on the three-dimensional model EFDC (Environmental Fluid Dynamics Code) and the HFR data using a Nudging or DI algorithm was considered the most appropriate for Galway Bay. The largest averaged Data Assimilation Skill Score (DASS) over the ≥6 h forecasting period from the best model NDA attained 26% and 31% for east–west and north–south surface velocity components respectively. Because of its ease of implementation and its accuracy, this data assimilation system can provide timely and useful information for various practical coastal hindcast and forecast operations.

1. Introduction

Globally, significant commercial and recreational activities occur along coasts; therefore, knowledge of coastal circulation patterns is of great importance. The timely information of surface flow patterns with high accuracy plays a key role in many coastal activities. Coastal process models are used extensively for hindcasting processes used in engineering projects, such as identifying the appropriate locations for sewage treatment plant sea outfall, or for investigating historical events such as accidental oil spillages. Hydrodynamic operational forecasting is becoming more common, and plays an important role in forecasting events such as storm surges, coastal flooding and for search and rescue missions. As electronic communication technology advances, installation of advanced remote sensing oceanic platforms based on radars or satellite mean that more and more near real-time ocean data such as surface currents and waves over a large domain and short observation windows are available [1,2,3]. A large number of radar observation systems have now been deployed in coastal areas around the world due to their relatively low cost and higher data density than satellite data sources.

Remote sensing data, especially surface currents, monitored by High-Frequency radar (HFR) system have been used to explore underlying dynamics of coastal waters, to validate models, to characterize coastal circulation patterns and to improve modelling performance [4,5]. Because synoptic, high spatial and temporal resolution surface flow fields are monitored in near real-time by the HFR system, application of HFR data in operational coastal models is becoming increasingly popular. Although coastal numerical models simulate states of parameters of interest, models do not always generate acceptable results; many of the difficulties with numerical models are due to poor definition of initial and boundary conditions. To enhance model performance a number of researchers have employed data assimilation techniques to combine available observations with models. Previous studies of various natural systems (such as hydrological systems) indicate that integration of accurate and efficient data assimilation algorithms into models can have significant benefits for model applications in both hindcasting and forecasting [6,7,8].

Data assimilation is a technique used to obtain a better estimate of parameter states by combining information of observations and models. It has been broadly applied in areas such as atmospherics, hydrology and oceanography [9,10,11,12]. Variational and sequential data assimilation algorithms are two popular approaches that have been used in oceanographic modelling. The various algorithms have varying degrees of complexity of implementation; HFR data of oceanic parameters can be assimilated into models using different data assimilation algorithms for different coastal zones. For example, Gopalakrishnan and Blumberg [6] assimilated HFR data into a regional model using a Nudging data assimilation algorithm in the New York Bight domain; Shulman and Paduan [13] combined HFR data into a model using an Optimal Interpolation data assimilation algorithm in the Monterey Bay area; Marmain et al. [14] used HFR data to optimize forcing in model using ensemble Kalman smoother in the northwestern Mediterranean Sea; Barth et al. [15] assimilated HFR data in a nested model of the West Florida Shelf using Ensemble Optimal Interpolation data assimilation algorithm; and Iermano et al. [16] applied four-dimensional variational data assimilation algorithm to assimilate HFR surface velocity fields, satellite sea surface height, sea surface temperature for the southern Tyrrhenian Sea. To date, little work has been carried out on inter-comparisons of a set of data assimilation algorithms and their underlying parameters under the same conditions for a given water body. Very little work has also been carried out on comparing the implementation of different data assimilation algorithms using HFR data in complex inshore waters, such as Galway Bay, which are heavily influenced by both wind and tidal dynamics [17,18,19].

With the goal of characterizing and predicting surface flows in Galway Bay area, a HFR system that generates hourly surface current flow fields of Galway Bay was used as a source for assimilating flow fields into a three-dimensional model EFDC (Environmental Fluid Dynamics Code) in this research [19,20,21], unlike toy models used artificial data such as Lorenz model [22]. One of the main objectives of this research was to assess the suitability of practical data assimilation algorithms (i.e., algorithms can be implemented relatively easily) for hindcasting and forecasting in shallow coastal waters subjected to alternate wetting and drying using data collected from a HFR system. To compare the efficacies of different data assimilation algorithms and to improve model simulations for Galway Bay, we tested four different algorithms: Direct Insertion (DI), Optimal Interpolation (OI), Nudging and indirect data assimilation via correcting model forcing into a hydrodynamic model and carrying out detailed inter-comparisons of model performances.

Models using each of the above data assimilation technique were implemented; specific model data assimilation parameters were based on sensitivity experiments of data assimilation cycle length times. Two statistical assessment methods, namely Data Assimilation Skill Score (DASS) and Averaged Kinetic Energy (AKE) were used to quantify the effects of data assimilation on modelling performance after employing the various data assimilation techniques. The best data assimilation model using HFR data was determined in terms of both quantitative and qualitative comparisons. Inter-comparison among the best four data assimilation models indicated that the hydrodynamic models using DI and Nudging improved modelling accuracy greater than hydrodynamic models using OI and indirect data assimilation.

The remainder of the paper is organized as follows: Section 2 introduces the methodology, including research domain, the numerical model, and a brief introduction about the HFR system deployed. Data assimilation algorithms are presented in Section 3, followed by results in Section 4. Discussion is presented in Section 5. Finally, conclusions are presented in Section 5.

2. Methodology

2.1. Research Domain

Galway Bay is located on the west coast of Ireland with the bay entrance opening onto the northeast Atlantic Ocean, as shown in Figure 1. The bay is semi-enclosed as it is partially shielded by the Aran Islands [23]. The water circulation through four sounds influences the physics, chemistry and biology of the waters of the bay [24,25]. Oceanic flows enter the bay mainly through the southern sounds; circulate within the bay before exiting through the North Sound [19]. The River Corrib drains approximately 70% of the catchment area around Galway Bay with a mean annual flow of 98.89 m³/s [26], discharging into the inner bay and generally extending seawards along the north coast [27]. The dynamics within the bay are mainly influenced by oceanic flows to the bay from the adjacent shelf and wind driven currents. Meteorological conditions in Galway Bay were mainly influenced by the Atlantic weather system [23]. This is typical of the prevailing weather conditions in Galway Bay. The prevailing wind direction in Galway Bay is from the southwestward. Waters in Galway Bay are well-mixed with no evidence of stratification [28]. Detailed review on hydrodynamic characteristics in this area, see Ren et al. [18] and O’Donncha et al. [19].

Galway Bay has one of the few open water marine renewable test sites in the world and so is of major local economic significance. The site has been structured to enable developers test quarter scale devices in the real environment. As part of this development, it is vital that there is a good understanding of the environmental conditions at the site. The authors have developed a detailed 3D hydrodynamic model to provide accurate hydrodynamic data at the site [18]. In this work, the computational domain of Galway Bay in the three-dimensional model EFDC as shown in Figure 1 is approximately 57 km in length and 36 km at its widest. The water depths of Galway Bay within the area covered by the HFR system are in the range of 10–40 m.

2.2. Numerical Model

The numerical model EFDC was used to simulate the hydrodynamic circulation of Galway Bay; the U.S. Environmental Protection Agency (EPA) developed the model at the Virginia Institute of Marine Science. EFDC comprises of four linked modules: hydrodynamics, water quality and eutrophication, sediment transport, and toxic chemical transport and fate. Only the hydrodynamic module was used during this research; this module solves the three-dimensional, vertically hydrostatic, free surface, turbulent averaged equations of motions for a variable density fluid. The hydrodynamic component of EFDC uses a semi-implicit, conservative finite volume solution scheme for the hydrostatic primitive equations with either two or three level time stepping [20,29,30]. The model uses a sigma vertical coordinate system and either rectangular grid, or curvilinear orthogonal horizontal coordinates. Over the past 25 years, the model has been applied to a variety of modelling studies of rivers, lakes, estuaries and coastal regions [31,32,33].

In this research, the EFDC model of Galway Bay (see Figure 1) in its barotropic mode was developed using a regular rectangular grid coordinate system; a 150 m horizontal spatial resolution was employed yielding 380 × 241 grid cells. A bathymetric model of Galway Bay was developed as part of the Integrated Mapping for the Sustainable Development of Ireland’s Marine Resource (INFOMAR) programme. The INFOMAR programme is an Irish national joint venture project between the Geological Survey of Ireland (GSI) (Dublin, Ireland) and the Irish Marine Institute (MI) (Galway, Ireland). The programme is the successor to the Irish National Seabed Survey (INSS) and concentrates on creating integrated mapping products related to the seabed of Irish territorial waters. The number of vertical layers in the EFDC model for Galway Bay was set to 10. Variable vertical layer thicknesses were used in the model with a thinner layer at the top and bottom of the water column and thicker layers in the middle, thereby ensuring that wind forcing was not overly-damped by tidal forcing. Detailed description on set up of the vertical layer structure for the Galway Bay is provided by Ren et al. [18]. The meteorological forcing parameters including wind, pressure, rain, solar radiation and relative humidity were obtained at one-minute intervals from the Informatics Research Unit for Sustainable Engineering (IRUSE) weather station located at the campus of National University of Ireland Galway, 4 km inshore from the northeast coast of Galway Bay. Because of its semi-enclosed nature, the dynamics in the Galway Bay area are a result of the oceanic flows to the bay from the adjacent shelf, inflows from the River Corrib and wind driven currents. Galway Bay usually experiences favourable conditions for marine resources. Tide and local, fetch-limited winds are the main driving forces for surface currents [34]. Records of the River Corrib inflows, which enter Galway Bay close to the north of point C1 in Figure 1, were obtained from the Irish Office of Public Works (OPW) (Trim, Ireland). Tidal water elevation time series generated from 2002 Oregon State University Tidal Inversion Software (OTIS) (Corvallis, OR, USA) were used to define the tidal forcing on the western and southern open boundaries in the model [35,36]. The tidal range during the analysis period of this work was approximately 5 m. The dominant tidal direction in Galway Bay area is in the east–west direction.

2.3. HFR System

A network of the Coastal Ocean Dynamics Application Radar (CODAR) SeaSonde HFR system had been deployed in Galway Bay since July 2011; they are located at Mutton Island Waste Water Treatment Plant and Spiddal Pier (see C1 and C2 in Figure 1). This HFR system remotely senses near-surface ocean currents and gravity surface waves. The system is based on the transmission of electromagnetic signals to the ocean and the receiving of a reflected signal by radio-wave backscatter, that is exactly equal to half of the transmitted signal wavelength [37,38]. Operating settings are: the transmission frequency is 25 MHz and the system bandwidth is 500 kHz. Summing radial data from two or more radars creates total surface vector fields at a 300 m horizontal resolution every 60 min. The coverage domain of the HFR system in the inner bay is highlighted in Figure 1. Data were filtered so that extreme values, i.e., greater than 100 cm/s, were not recorded, as it is known from historic observations 100 cm/s is an upper limit on current velocities in Galway Bay. Along with the surface velocity components, the radar systems also determine the standard deviations of both surface velocity components. The averaged standard deviations are 1.01 cm/s and 0.84 cm/s for east–west and north–south surface velocity component during Julian Day 182–201 in 2013, respectively. The range of HFR surface velocity components is approximately 0–100 cm/s. The radar data standard deviation is 1.1% on average. This indicates that HFR system can stably monitor surface currents within acceptable limits.

The accuracy of the HFR system along the baseline is conceivably lower because of geometrical configuration issues, known as geometric dilution of precision (GDOP) [39]. GDOP is a coefficient of uncertainty that characterizes the effect of the geometry of the coupled radar system on the measurement and position determination errors. Near the baseline, total vectors are not accurate as the radial velocities are nearly parallel [19]. O’Donncha et al. [19] had computed Pearson’s coefficient between the HFR data and the Galway Bay model “Free Run” without employing data assimilation (see Figure 2). Model skill scores proposed by Willmott [40] were also computed to better quantify localized discrepancies allowing for ambient flow speeds (see Figure 3), and in particular lower flows in the baseline region along the Northern shore. Both metrics are presented decomposed into the zonal and meridional component of flow.

Figure 2 and Figure 3 show that the meridional component of flow is distorted most by GDOP along the baseline, while the zonal component is accurately resolved, apart from a very small domain along the shoreline because of a slight rotation of the baseline from east–west. Analysing these results, Figure 2b and Figure 3b suggest that GDOP does not have a significant impact on results. There is no evidence of a decrease in agreement between the meridional datasets along this narrow domain [19]. Measurements of surface currents in Galway Bay have been validated with independent observations from an Acoustic Doppler Current Profile (ADCP), see details by O’Donncha et al. [19] and Ren et al. [41]. In this work, HFR surface currents were assimilated into an EFDC model of Galway Bay.

The radar system area coverage includes extensive commercial uses such as fishery, tourism, aquaculture and marine energy. Shore-based radars measure high resolution spatial and temporal evolution of the two-dimensional surface current field over relatively large coastal domains [42]. However, the domain covered by the HFR system changes over time due to variable roughness of ocean surface. Because the CODAR system can provide continuous dataset over years once it is installed and well maintained, the understanding, extraction and application of these data are new challenges to explore the dynamical process of surface currents. The synoptic, high resolution, near real-time view of surface flow fields makes it possible to study local dynamics and transport processes over a large domain for years. Thus HFR systems have a number of useful real world applications: HFR data are used in search and rescue operations by US Navy Coast Guard [43]; HFR data are used to improve oil spill modelling performance [44,45,46]; HFR data are assimilated into models to enhance modelling accuracy [14,47,48,49]. In this work, HFR surface currents were assimilated into models to explore potential for improving modelling performance, and to assess the accuracy of several potential data assimilation techniques.

3. Data Assimilation Algorithms

Two generic types of data assimilation algorithms are widely used in oceanic forecasting systems: firstly, sequential data assimilation algorithms such as OI, DI and Nudging, in which observations are linearly combined into model background states to produce the optimized analysis states; and, secondly, variational data assimilation algorithms such as three-dimensional (3D) and four-dimensional (4D) variational algorithms, in which analysis states are obtained through a cost function. The cost functions sum two distances: one is the distance between the analysis states and observation states, and the other is the distance between the analysis states and the model background states. Only sequential data assimilation techniques were used in this research. A complete data assimilation system is generally comprised of three components: (a) a numerical model, which describes the dynamic process in a mathematic way; (b) a batch of observations of variables of interest; and (c) a data assimilation algorithm which blends the measurement states with the model background states [50,51].

In this work, the 3D numerical coastal EFDC model was adapted to integrate data assimilation procedures; surface currents measured by the CODAR system were assimilated into the model to enhance its performance. HFR data were bilinearly interpolated on to the higher resolution model grid before data assimilation. As outlined above, suitability of sequential data assimilation or use in the Galway Bay EFDC forecasting/hindcasting model was assessed.

The assessment can be considered in two stages: firstly, the assimilation parameters of each data assimilation technique were examined through sensitivity analysis; and secondly, data assimilation cycle lengths were examined and determined using temporally interpolated HFR data. Detailed description of each data assimilation algorithm is presented in the following sections.

For simplicity, the following notation is used in this paper: (a) DI, data assimilation using the Direct Insertion algorithm [52]; (b) OI, data assimilation using the Optimal Interpolation algorithm [53]; (c) NDA, data assimilation using Nudging algorithm [54]; (d) IDA, data assimilation using indirect data assimilation algorithm through correcting wind stress [41]; and (e) FR, a “Free Run”, a benchmark model without data assimilation implemented. All the initial and boundary conditions are the same in all the above models. The simulation period was from Julian Day 211 to 230 in 2013; this period was chosen as the measured HFR data had high spatial density during this period relative to other time periods. This period is subdivided as: (i) spin-up period, Julian Day 211–220; (ii) data assimilation period, Julian Day 220–228 01:00; and (iii) forecasting period, after Julian Day 228 01:00. Introduction of each used data assimilation technique is presented in the following sections.

3.1. OI DA Algorithm

The analysis equation of the OI data assimilation algorithm can be expressed as [55]:

$$\widehat{{\mathrm{x}}^{\mathrm{a}}}=\widehat{{\mathrm{x}}^{\mathrm{b}}}+\mathrm{K}\left(\widehat{{\mathrm{y}}^{\mathrm{o}}}-\mathrm{H}\widehat{{\mathrm{x}}^{\mathrm{b}}}\right)$$

(1)

where $\widehat{{\mathrm{x}}^{\mathrm{a}}}$ is the analysis state; $\widehat{{\mathrm{x}}^{\mathrm{b}}}$ is the model background state; $\mathrm{K}$is the Kalman gain; $\widehat{ {\mathrm{y}}^{\mathrm{o}}}$is the measurement state; and $\mathrm{H}$is the observation projector.

The Kalman gain in OI is derived from minimum analysis error covariance [56]:

$$\mathrm{K}=\frac{{\mathrm{P}}^{\mathrm{b}}{\mathrm{H}}^{\mathrm{T}}}{{\mathrm{HP}}^{\mathrm{b}}{\mathrm{H}}^{\mathrm{T}}+\mathrm{R}}={\mathrm{P}}^{\mathrm{b}}{\mathrm{H}}^{\mathrm{T}}{\left({\mathrm{HP}}^{\mathrm{b}}{\mathrm{H}}^{\mathrm{T}}+\mathrm{R}\right)}^{-1}$$

(2)

$$\mathrm{R}=\mathrm{E}({\mathsf{\epsilon }\mathsf{\epsilon }}^{\mathrm{T}})=\left[\begin{array}{ccc}\mathrm{var}\left({\mathsf{\epsilon }}_1\right)& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & \mathrm{var}\left({\mathsf{\epsilon }}_{\mathrm{n}}\right)\end{array}\right]$$

(3)

$${\mathrm{P}}^{\mathrm{b}}=\mathrm{E}\left(\widehat{{\mathsf{\eta }}^{\mathrm{b}}}{\widehat{{\mathsf{\eta }}^{\mathrm{b}}}}^{\mathrm{T}}\right)=\left[\begin{array}{ccc}\mathrm{var}\left({\mathsf{\eta }}_1^{\mathrm{b}}\right)& \cdots & \mathrm{cov}\left({\mathsf{\eta }}_1^{\mathrm{b}},{\mathsf{\eta }}_{\mathrm{n}}^{\mathrm{b}}\right)\\ ⋮& \ddots & ⋮\\ \mathrm{cov}\left({\mathsf{\eta }}_{\mathrm{n}}^{\mathrm{b}},{\mathsf{\eta }}_1^{\mathrm{b}}\right)& \cdots & \mathrm{var}\left({\mathsf{\eta }}_{\mathrm{n}}^{\mathrm{b}}\right)\end{array}\right]$$

(4)

$$\mathrm{cov}\left({\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}},{\mathsf{\eta }}_{\mathrm{j}}^{\mathrm{b}}\right)=\sqrt{\mathrm{E}({\left({\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}}\right)}^2)-{\left[\mathrm{E}({\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}})\right]}^2}\times \mathsf{\rho }({\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}},{\mathsf{\eta }}_{\mathrm{j}}^{\mathrm{b}})\times \sqrt{\mathrm{E}({({\mathsf{\eta }}_{\mathrm{j}}^{\mathrm{b}})}^2)-{\left[\mathrm{E}({\mathsf{\eta }}_{\mathrm{j}}^{\mathrm{b}})\right]}^2}$$

(5)

$$\mathsf{\rho }\left({\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}},{\mathsf{\eta }}_{\mathrm{j}}^{\mathrm{b}}\right)=\exp (-\frac{{\mathrm{r}}^2}{{\mathrm{L}}^2})$$

(6)

where ${\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}}$ is the model background error at the ith grid point, ${\mathsf{\eta }}_{\mathrm{i}}^{\mathrm{b}}={\mathrm{x}}_{\mathrm{i}}^{\mathrm{b}}-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}}$; ${\mathsf{\epsilon }}_{\mathrm{i}}$ is the measurement error at the ith grid point, ${\mathsf{\epsilon }}_{\mathrm{i}}={\mathrm{y}}_{\mathrm{i}}^0-{\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}}$; ${\mathrm{x}}_{\mathrm{i}}^{\mathrm{t}}$ is the true state at the ith grid point; E(.) is the expectation; $\mathrm{cov}(.)$ is covariance; L is the horizontal correlation length; r is the distance between ith and jth grid points; $\mathrm{T}$ is the matrix transpose; ${\mathrm{P}}^{\mathrm{b}}$ is the model background error covariance; and $\mathrm{R}$ is the observation error covariance.

In the OI data assimilation algorithm, the key aspect is the definition of the Kalman gain $\mathrm{K}$. In this work, the measurement error covariance was the error covariance of the HFR data. Since the Galway Bay HFR radar system produces standard deviations (SD) for the surface velocity components at each measurement step, the measurement error covariance matrix $\mathrm{R}$ was determined from the measurement error variances at each data grid point. Based on Martin et al. [57] using the pairs of observations minus model background states to form model background error statistics for assimilating ocean parameters such as sea surface temperature and salinity, similar formulations from HFR observations minus model background states were used to define model background errors in this work.

3.2. DI DA Algorithm

The DI algorithm is essentially a simplified version of the OI algorithm. The analysis states are replaced by the measurement states directly. In other words, no consideration of model errors and measurement errors was required leading to the Kalman gain being an identity matrix. The analysis equation is simplified to [52,58]:

$$\widehat{{\mathrm{x}}^{\mathrm{a}}}=\widehat{{\mathrm{y}}^{\mathrm{o}}}$$

(7)

At the model grid points where available observations exist, the model background states are directly replaced by the measurement states. A major assumption that measurement states are viewed as “perfect” is implicit in the DI algorithm. Since the EFDC model has not been used to assimilate HFR surface currents, application of the DI data assimilation algorithm in this work is used to examine robustness and potential of using EFDC model to develop an operational data assimilation forecasting system.

3.3. Nudging DA Algorithm

The nudging data assimilation algorithm is quite similar to the OI algorithm. It is developed by adding the weighted differences between the measurement states and model background states to the model background states. However, the definition of the weighting factor as specified as the Kalman gain in OI was different from that used in the Nudging algorithm. The Kalman gain, K, of Equation (1) was replaced by a nudging parameter $\mathsf{\lambda }$. In coastal data assimilation models, Fan et al. [59], Lin et al. [60] and Gopalakrishnan [61] defined the Nudging data assimilation analysis equation as:

$$\frac{\partial \mathrm{U}}{\partial \mathrm{t}}=\left(\mathrm{physics}\right)+\left[\frac{1}{{\mathrm{t}}_{\mathrm{a}}}\right]{\ast \mathrm{e}}^{\left[\frac{-{\mathrm{r}}_{\mathrm{nudging}}^2}{{\mathrm{d}}_{\mathrm{n}}^2}\right]}{\ast \mathrm{e}}^{\left[\frac{-\left(\mathrm{t}-{\mathrm{t}}_0\right)}{{\mathrm{t}}_{\mathrm{d}}}\right]}{\ast \mathrm{e}}^{\left[\frac{\mathrm{z}}{{\mathrm{z}}_{\mathrm{d}}}\right]}\ast \left({\mathrm{U}}^0-\mathrm{U}\right)$$

(8)

where U is the model background states; ${\mathrm{U}}^0$ is the observation states from radars; (physics) denotes the physical process mathematically described in the numerical model; ${\mathrm{r}}_{\mathrm{nudging}}$ is the distance between model grid point and the observation location; $\left(\mathrm{t}-{\mathrm{t}}_0\right)$is the difference between assimilation and observation time; ${\mathrm{t}}_{\mathrm{a}}$is the assimilation timescale, which determines the strength of the nudging parameter; ${\mathrm{t}}_{\mathrm{d}}$is the damping time-scale for the nudging term; ${\mathrm{d}}_{\mathrm{n}}$ is the length-scale of nudging term; ${\mathrm{e}}^{\left[\frac{\mathrm{z}}{{\mathrm{z}}_{\mathrm{d}}}\right]}$is an exponential decay parameter, which controls the depth of influence of the nudging parameter; and ${\mathrm{z}}_{\mathrm{d}}$ is the depth of influence.

A value of ${\mathrm{e}}^{\left[\frac{-{\mathrm{r}}_{\mathrm{nudging}}^2}{{\mathrm{d}}_{\mathrm{n}}^2}\right]}$equal to 1 indicates that the model grid point and the observational location are the same, and ${\mathrm{e}}^{\left[\frac{-\left(\mathrm{t}-{\mathrm{t}}_0\right)}{{\mathrm{t}}_{\mathrm{d}}}\right]}=1$means that there is no time-lag $\left(\mathrm{t}-{\mathrm{t}}_0=0\right)$ between model and measurements. If the model grid points are close to observational data points and there is no time-lag at data assimilation step for measured data, the nudging parameter can be simplified in terms of the two conditions [61]:

$$\mathsf{\lambda }=\left[\frac{1}{{\mathrm{t}}_{\mathrm{a}}}\right]{\ast \mathrm{e}}^{\left[\frac{\mathrm{z}}{{\mathrm{z}}_{\mathrm{d}}}\right]}$$

(9)

The nudging parameter $\mathsf{\lambda }$ is zero at grid points, where HFR observations are not available. Since interpolated HFR data accord with model grid and HFR data were temporally interpolated on to each analysis model computational time step, the same nudging as described in Equation (9) was used in this work. Thus, by substituting Equation (9) into Equation (1), the analysis equation for the Nudging data assimilation algorithm is expressed as:

$$\widehat{{\mathrm{x}}^{\mathrm{a}}}=\widehat{{\mathrm{x}}^{\mathrm{b}}}+\left[\frac{1}{{\mathrm{t}}_{\mathrm{a}}}\right]{\ast \mathrm{e}}^{\left[\frac{\mathrm{z}}{{\mathrm{z}}_{\mathrm{d}}}\right]}\left(\widehat{{\mathrm{y}}^{\mathrm{o}}}-\mathrm{H}\widehat{{\mathrm{x}}^{\mathrm{b}}}\right)$$

(10)

Equation (10) was implemented in the Galway Bay model by Ren and Hartnett [54].

3.4. Indirect DA Algorithm

Because patterns of surface currents in the research domain are strongly affected by wind forces, it is possible that model errors result from inaccurate wind stress specification. To take advantages of the HFR data and to improve simulation of strongly wind-induced surface currents, HFR data were used to correct wind stresses in the models. In this approach, firstly, differences of surface currents between model results and the HFR data were computed; then the differences at each measurement step were adopted to update wind stress using the algorithm proposed by Lewis et al. [62]. Detailed description of this approach for Galway Bay domain is presented by Ren et al. [41].

4. Results

The best assimilation parameters in each data assimilation algorithm were selected based on sensitivity tests generating closer patterns of surface currents to the HFR data during hindcasting period. For OI system, the best horizontal correlation length was found to be 50 km [58]. For Nudging data assimilation system, the best depth of influence and time scale were found to be 4.0 m and 1800 s, respectively [54]. For indirect data assimilation via correcting wind stress, the best assimilation coefficient is −1 [41]. The Authors found that hourly assimilation of the HFR data did not have significantly positive impacts on forecasting when using these data assimilation algorithms. To explore the potentials of for enhancing model forecasting capability, sensitivity tests on data assimilation cycle lengths were examined in each algorithm separately. Data assimilation cycle lengths less than the HFR data observation interval 1 h (15 min, 5 min, 1 min and model computational time step) were selected and examined. Other data assimilation cycle lengths less than 1 h were examined as well, but the differences were not significant. The main reasons for assimilating data over shorter temporal intervals into models are presented by Ren and Hartnett [54]. The best data assimilation cycle length was found to be at each model computational time step for DI, OI and Nudging; the best data assimilation cycle length was one minute for indirect data assimilation via correcting wind stress. For details of the procedure using with IDA and OI data assimilation algorithm, see Ren et al. [41] and Ren et al. [58], respectively. The same procedure was used to determine data assimilation cycle length for assimilation models using the Nudging algorithm by Ren and Hartnett [54].

In this work, inter-comparisons among the various data assimilation models were undertaken to assess their accuracy during both hindcasting and forecasting, and to select the most suitable data assimilation model for Galway Bay based on this comprehensive assessment. Detailed comparisons are presented in the following sections. Since in the Direct Insertion data assimilation algorithm model background states were replaced directly with HFR data during assimilation period, analysis states were the same as HFR data. Direct Insertion was primarily applied to examine the robustness of the EFDC for combining radar data into models. Surface vector fields from DI models during hindcasting period are thus not shown in Section 4.1.

4.1. Hindcasting of Surface Currents

Because the Authors focused on simulating patterns of surface currents for Galway Bay, vector fields of surface currents at a typical time (Julian Day 224 01:00) from the best data assimilation model using each of the algorithms along with HFR data are presented in Figure 4.

Figure 4 shows that assimilation models, especially model OI and NDA, had better agreement with the HFR data (Figure 4e) than model FR. For completeness, it is instructive to compare how the various models behave during hindcast mode as well as forecast mode. To clearly show the improvements from the best assimilation models, results from model DI were not shown in Figure 4. Southeast flow trends near the northern corner of the simulation domain were obtained in models OI and NDA in comparison with HFR data, but models IDA and FR did not capture them. Moreover, model OI only captured similar surface flow trends as the HFR data for water that flows away from the coast near radar station C1. The Root-Mean-Square-Error (RMSE) values between models and radar data at Julian Day 224 01:00 are shown in

Table 1. RMSE between models and radar data (Julian Day 224 01:00).

Model	RMSE (u) (cm/s)	RMSE (v) (cm/s)
FR	6.46	6.50
OI	1.48	0.96
NDA	5.65	4.61
IDA	6.18	6.20

.

Improvements in RMSE values between results from model OI and HFR data for surface east–west and north–south velocity component in comparison with model FR are 77% and 85%, respectively, at Julian Day 224 01:00, while improvements in RMSE values for model NDA are 13% and 29% for surface east–west and north–south velocity component, separately. Improvements in RMSE values for surface east–west and north–south velocity components in model IDA are not as significant as the other models at 4% and 5%, respectively. The above results indicate that all data assimilation models improved modelling performance during hindcasting period. The OI algorithm forces model background states closer to HFR observations.

4.2. Forecasting of Surface Currents

Accurate forecasting of sea states is of great importance for many coastal operations such as search and rescue and shipping guidance. To assess the patterns of surface currents during the forecasting period, vector fields of surface current at two representative forecasting times Julian Day 228 02:00 and 03:00 from the models and the HFR data are shown in Figure 5 and Figure 6, respectively.

Figure 5 shows that the assimilation models were closer to the HFR data (Figure 5f) than the model without data assimilation (FR) at Julian Day 228 02:00. Improvements in RMSE between results from assimilation models and HFR data in model NDA and DI are significant when compared with model FR. Model NDA improved by 10% and 37%, respectively, for surface east–west and north–south velocity components; and model DI by 12% and 37%, respectively, for surface east–west and north–south velocity components. However, the improvement in the RMSE of the surface east–west velocity component in model OI is more significant by 17%, whereas the RMSE of the surface north–south velocity component is improved by 10%. The improvement in the RMSE for the surface north–south velocity component in model IDA is 22%; however, data assimilation caused degradation of the simulation of the east–west surface velocity components in model IDA. This is likely to have resulted due to the application of a constant wind stress during the assimilation of HFR data over the domain; wind stresses actually varied in space over the domain.

Figure 6 shows that the assimilation models were in relatively good agreement with the HFR data (Figure 6f) in comparison with model FR at Julian Day 228 03:00. Improvements in RMSE between results from assimilation models NDA and DI, and HFR data are comparable. Model NDA RMSE values improved by 10% and 34%, respectively, for surface east–west and north–south velocity component; and model DI RMSE values improved by 11% and 33%, respectively, for surface east–west and north–south velocity component in comparison with model FR. However, the improvement in the RMSE of the east–west surface velocity component in model OI is more significant by 20% than the north–south surface velocity component by 6%. The improvement in RMSE for north–south surface velocity component in model IDA is 25%; however, data assimilation caused degradation of the simulation of east–west surface velocity component in model IDA. This likely resulted again from using a constant wind stress update over the domain.

To show the temporal changes of surface velocity components during forecasting, time series of surface velocity components at point A (see Figure 1) during the ≥6 h forecasting (Julian Day 228 02:00–07:00) are shown in Figure 7 and Figure 8.

Figure 7 and Figure 8 show that forecasting of both surface velocity components was improved in each of the best data assimilation models over the ≥6 h forecasting compared with the model FR. For the east–west velocity components, results from data assimilation models had closer trends to HFR data than results from model FR. Additionally, the data assimilation models were comparable with each other, however, model IDA outperformed other data assimilations at certain times. For the north–south velocity components, the data assimilation models were closer to HFR data than model FR. Results from models DI, IDA and NDA were closer to HFR data than model OI over the ≥6 h forecasting period.

To illustrate the degree of improvement of the data assimilation models over model FR, spatial variations in RMSE between radar data and model results were calculated using Equations (11) and (12) and are presented in Figure 9. RMSE of surface velocity east–west components (u) were firstly calculated at each assimilation point using Equation (11) during the ≥6 h forecasting period (Julian Day 228 02:00–07:00). Equation (11) was also used to compute RMSE values for north–south surface velocity components (v). The total $\mathrm{RMSE}\left(\mathrm{u},\mathrm{v}\right)$ at each assimilation point was then computed using Equation (12). The distribution of $\mathrm{RMSE}\left(\mathrm{u},\mathrm{v}\right)$ across the domain was used to assess the degree of agreement between modelled results and the HFR data.

$$\mathrm{RMSE}\left(\mathrm{u}\right)=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{j}}}{\left({\mathrm{u}}_{\mathrm{CODAR}}\left(\mathrm{i}\right)-{\mathrm{u}}_{\mathrm{model}}\left(\mathrm{i}\right)\right)}^2}{{\mathrm{N}}_{\mathrm{j}}}}$$

(11)

$$\mathrm{RMSE}\left(\mathrm{u},\mathrm{v}\right)=\sqrt{\mathrm{RMSE}{\left(\mathrm{u}\right)}^2+\mathrm{RMSE}{\left(\mathrm{v}\right)}^2}$$

(12)

where ${\mathrm{u}}_{\mathrm{CODAR}}$is the CODAR HF radar measured surface velocity east–west component (cm/s); ${\mathrm{u}}_{\mathrm{model}}$ is the modelled surface velocity east–west component (cm/s); $\mathrm{RMSE}(\mathrm{u}$), $\mathrm{RMSE}\left(\mathrm{v}\right)$ is the RMSE of east–west velocity and north–south component at time step $\mathrm{j}$ respectively (cm/s); ${\mathrm{N}}_{\mathrm{j}}$ is the number of calculation time steps; and $\mathrm{RMSE}\left(\mathrm{u},\mathrm{v}\right)$ is root squared RMSE value.

Figure 9 shows distributions of RMSE values over the data assimilation domain during ≥6 h forecasting period. RMSE values were generally smaller for model DI (Figure 9b) and model NDA (Figure 9d) than model FR (Figure 9a), model OI (Figure 9c) and model IDA (Figure 9e), especially in the area with high HFR coverage density. This indicates that, when assimilation models are rich in HFR data, these models improve forecasts, especially when using the DI or NDA algorithms.

4.3. Data Assimilation Skill Score

To quantify improvements in accuracy of modelled surface currents achieved by using different data assimilation techniques, a data assimilation skill score (DASS), which is based on the Mean-Square-Error (MSE), is calculated over time in the domain covered by the radar system. The DASS can be expressed as [63]:

$$\mathrm{DASS}=1-\frac{\mathrm{MSE}\left({\mathrm{v}}^{\mathrm{CODAR}},{\mathrm{v}}^{\mathrm{DA}}\right)}{\mathrm{MSE}\left({\mathrm{v}}^{\mathrm{CODAR}},{\mathrm{v}}^0\right)}$$

(13)

where ${\mathrm{v}}^{\mathrm{DA}}$ is the output (east–west or north–south velocity component) from the model with data assimilation; ${\mathrm{v}}^0$ is the output from the model FR; and ${\mathrm{v}}^{\mathrm{CODAR}}$ is the data from the HF radar system.

If DASS is greater than zero, then the data assimilation model improves the accuracy of the forecast over model FR. If DASS is less than zero, then the data assimilation contaminates the basic dynamic process in the model resulting in deterioration in model accuracy. Averaged DASS of surface velocity components from the best data assimilation model using each data assimilation technique were separately calculated over the ≥6 h forecasting (Julian Day 228 02:00–07:00). DASS for both surface velocity components over time indicated that assimilation models improved the forecasting with positive values except for the east–west velocity component in model IDA. The largest averaged DASS (east–west component: +26%, north–south component: +31%) over the ≥6 h forecasting was from the best model NDA. The magnitudes of DASS obtained here were greater than the results (east–west component of M2 mooring: +6%, north–south component of M2 mooring: +14%; east–west component of M3 mooring: +26%, north–south component of M3 mooring: +10%) in studies by Gopalakrishnan [61]. Significantly, Gopalakrishnan [61] developed DASS at a few discrete mooring points only; whereas in this research DASS was a domain averaged value, indicating substantial domain-wide improvements due to data assimilation.

To further compare DASS values for surface velocity components between the various data assimilation models, DASS values were averaged for + 0–3 h and + 4–6 h forecasting period, respectively, and shown in

Table 2. Averaged DASS values for surface velocity components during forecasting period.

Model	East–West Component (u)		North–South Component (v)
	+ 0–3 h	+ 4–6 h	+ 0–3 h	+ 4–6 h
DI	0.20	2.30	0.47	0.17
OI	0.28	0.15	0.09	−0.02
NDA	0.20	0.33	0.47	0.18
IDA	−0.13	0.15	0.34	0.12

.

Averaged DASS magnitudes in Table 2 shows that DASS decreased with time more significantly for north–south (v) surface velocity components than east–west surface velocity components (u). This is probably because tidal forcing alone is more dominant in the east–west direction. Forecasts of surface flows in the north–south direction, which are more strongly affected by winds, are shown to be improved through data assimilation. DASS values of models DI and NDA had bigger variation from + 0–3 h to + 4–6 h than models OI and IDA for north–south surface velocity components. Averaged DASS values of north–south surface velocity components in model OI were close to zero. This indicates that the effects of assimilating HFR data into models were not significant. However, the decreasing trend of averaged DASS values in model OI was the most significant for east–west surface velocity components. Models DI and NDA had an increasing trend for east–west surface velocity component from + 0–3 h to + 4–6 h.

4.4. Average Kinetic Energy

To further quantitatively assess the improvements of data assimilation in the hydrodynamic models using radar surface currents, spatially Averaged Kinetic Energy (AKE) $\overline{{\mathrm{E}}_{\mathrm{k}}}$ in the domain covered by the radar system for both model results and observations was calculated using the following [47,64]:

$$\overline{{\mathrm{E}}_{\mathrm{k}}}=\frac{1}{\mathrm{A}}\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}{\mathrm{U}}^2\mathrm{ds}$$

(14)

where $\mathrm{A}$is the domain area covered by the radar; and $\mathrm{U}$is the horizontal total current speed.

Water density was assumed constant and so omitted from Equation (14). Firstly, the AKE from models and observations were computed at each analysis time step. Hourly data over a 48-h forecasting period (Julian Day 228 01:00–Julian Day 230 0:00, 2013) were used for this analysis, so there are 48 AKE values used from each dataset to ensure a meaningful analysis. Then correlations of AKE between model results and the radar observations were computed. To show the influence of data assimilation over a reasonably long forecasting period, the correlation of AKE during a two-day forecasting is presented in

Table 3. Correlation of AKE of surface vector speeds.

Model	FR vs. CODAR	DI vs. CODAR	OI vs. CODAR	NDA vs. CODAR	IDA vs. CODAR
Cor(AKE)	0.52	0.59	0.54	0.59	0.54

.

The correlation coefficients of AKE between the HFR data and the models results in the above table indicate that stronger correlation appeared in all assimilation models than in model FR. The models NDA and DI had better correlation with the HFR data than models IDA and OI. The improvement in AKE correlation in comparison against model FR was 13% for the models NDA and DI. The improvement in AKE correlation in models OI and IDA was smaller by 4%. Improvement in AKE values indicates that employment of data assimilation using HFR data enhanced model forecasting over the domain for a long period; models NDA and DI outperformed models OI and IDA.

4.5. Computational Cost

The computational times for the data assimilation models using DI, Nudging and indirect data assimilation via correcting wind stress were very similar under the same conditions on the Irish Centre for High-End Computing (ICHEC) platform; it took around 28 h for full simulations of hindcasting and forecasting. While the best model using the OI data assimilation algorithm took around 85 h for the same simulation period. Although OI data assimilation algorithm has a more detailed theoretical background, the performance and computational cost of the best OI data assimilation model was not as good as the simpler data assimilation models using Nudging and DI algorithm. More efficient computation platforms or computational approaches need to be considered when using OI data assimilation algorithm, especially when frequent data assimilation updating, such as at each model computational time step, is required.

5. Discussion

Four types of sequential data assimilation schemes—DI, OI, Nudging and indirect data assimilation via correcting model wind stress—were used to assimilate HFR surface current flow fields into a three-dimensional EFDC model for Galway Bay. No other researchers have carried out such as extensive assessment and inter-comparison of data assimilation techniques of a coastal hydrodynamic model. Model DI produced the same results as the HFR data because the model background states were directly replaced by the HFR data during hindcasting period. Each of the other three data assimilation schemes improved model hindcasting performance in terms of the values of RMSE. The improvements obtained during this research were much greater than the improvements obtained by Ragnoli et al. [65] (east–west component: +1%, north–south component: +0.02%); their values were obtained by averaging their RMSE values at ten points during the entire simulation period when using an OI algorithm. RMSE values by Ragnoli et al. [65] were not calculated during the forecasting period, but during the whole simulation period including spin-up period, assimilation period and forecasting period. Each of the best data assimilation models in this research improved the model forecasting based on RMSE values. The best Nudging and DI data assimilation model, which assimilated the HFR data at each model computational time step, generated the largest improvements compared with other best assimilation models (OI and IDA). The best Nudging data assimilation model was reasonably accurate and quite efficient for the research domain. Values of RMSE throughout the domain with high HFR coverage density during ≥6 h forecasting period were smaller in general when using model DI and model NDA than the RMSE values obtained using the other assimilation models and model FR.

Data Assimilation Skill Score for both surface velocity components over time indicated that assimilation models improved forecasting with positive DASS values except for the east–west velocity component in model IDA. The largest averaged DASS values (east–west component: +26%, north–south component: +31%) over the ≥6 h forecasting were from the best model NDA. The magnitudes of DASS were greater than results obtained in the studies by Gopalakrishnan [61]. Significantly, Gopalakrishnan [61] only developed DASS at a few points, their mooring, whereas in this research DASS was a domain averaged value, indicating substantial domain-wide improvements due to data assimilation. The four best assimilation models (DI, OI, NDA and IDA) had stronger correlations of Averaged Kinetic Energy with the HFR data over + 48 forecasting compared with model FR. The largest improvement (13%) of AKE was from the best model NDA and DI. The improvement of AKE correlation was comparable to an improvement of 16% during a ≥2 h forecast in that study by Breivik and Sætra [47]. The computational times for the best models using DI, Nudging and indirect data assimilation via correcting wind stress were very similar under the same conditions on the Irish Centre for High-End Computing (ICHEC) platform; it took around 28 h for full simulations, while the best model using the OI data assimilation algorithm took around 85 h for a similar simulation.

6. Conclusions

The main conclusions arising from the research presented in this paper are as follows.

(1)

The three-dimensional EFDC model is robust enough to frequently combine measurements from the coastal radar system into routine data assimilation algorithms in complex inshore waters strongly influenced by both tides and wind dynamics. Forecasting improvements: Each of the best data assimilation models improved the model forecasting based on the values of RMSE. The best Nudging and DI data assimilation model, which assimilated the HFR data at each model computational time step, generated the largest improvements compared with other best assimilation models (OI and IDA). The best Nudging data assimilation model was reasonably accurate and quite efficient for the research domain. Values of RMSE throughout the domain with high HFR coverage density during ≥6 h forecasting period were smaller in general when using model DI and model NDA than the RMSE values obtained using the other assimilation models and model FR.

(2)

All data assimilation models improved modelling performance during a hindcasting period. The OI algorithm was found to force model background states closer to HFR observations.

(3)

Application of those data assimilation algorithms using radar data improved forecasts of north–south surface velocity components to a greater degree than east–west surface velocity components, while model IDA degraded simulation of east–west surface velocity components. This likely resulted from using a constant wind stress over the entire computational domain.

(4)

Forecasting of east–west surface velocity components from all data assimilation models were closer to HFR data than model FR. Results from data assimilation models in general were comparable with each other. For the north–south surface velocity components, results from models DI, IDA and NDA were closer to HFR data than results from model OI over a ≥6 h forecasting period.

(5)

Distributions of RMSE values over the data assimilation domain during ≥6 h forecasting period indicate that when assimilation models are rich in HFR data these models improve forecasts, especially when using the DI or NDA algorithms.

(6)

DASS values decreased with time more significantly for north–south surface velocity components than east–west surface velocity components. This is probably because tidal forcing alone is more dominant in the east–west direction. Moreover, improvement in AKE values indicates that employment of data assimilation using HFR data enhanced model forecasting over the domain for a long period; models NDA and DI were shown to outperform models OI and IDA.

In summary, sequential data assimilation algorithm is a powerful tool to enhance model forecasting of surface currents using the HFR data. Determination of appropriate assimilation parameters and data assimilation cycle lengths is crucial for producing meaningful forecasts. The NDA model using Nudging data assimilation algorithm and DI model using Direct Insertion data assimilation algorithm, were seen to produce the best forecasting results for Galway Bay area from accuracy and computational efficiency perspectives. Finally, it is important to remember that data assimilation parameters must be tuned for particular study domains.