Hindcasting and Forecasting of Surface Flow Fields through Assimilating High Frequency Remotely Sensing Radar Data

In order to improve the forecasting ability of numerical models, a sequential data assimilation scheme, nudging, was applied to blend remotely sensing high-frequency (HF) radar surface currents with results from a three-dimensional numerical, EFDC (Environmental Fluid Dynamics Code) model. For the first time, this research presents the most appropriate nudging parameters, which were determined from sensitivity experiments. To examine the influence of data assimilation cycle lengths on forecasts and to extend forecasting improvements, the duration of data assimilation cycles was studied through assimilating linearly interpolated temporal radar data. Data assimilation nudging parameters have not been previously analyzed. Assimilation of HF radar measurements at each model computational timestep outperformed those assimilation models using longer data assimilation cycle lengths; root-mean-square error (RMSE) values of both surface velocity components during a 12 h model forecasting period indicated that surface flow fields were significantly improved when implementing nudging assimilation at each model computational timestep. The Data Assimilation Skill Score (DASS) technique was used to quantitatively evaluate forecast improvements. The averaged values of DASS over the data assimilation domain were 26% and 33% for east–west and north–south velocity components, respectively, over the half-day forecasting period. Correlation of Averaged Kinetic Energy (AKE) was improved by more than 10% in the best data assimilation model. Time series of velocity components and surface flow fields were presented to illustrate the improvement resulting from data assimilation application over time.


Introduction
Accurate forecasting of surface currents in coastal areas is of great importance for operations such as search and rescue, fishing, and pollution monitoring. There are generally three main approaches to determine dynamic processes of coastal waters: observations, scaled physical modelling, and numerical modelling. Each of these provides meaningful information to analyze dynamic processes of relevant phenomena. However, each approach has its own constraints and shortcomings. Observations based on remote sensing are a desirable way to obtain data of real states, but they are usually expensive. Moreover, observation biases from observers or measurement tools can reduce the reliability of recorded data. Scaled physical experimentation is a useful way to investigate variations of dynamic processes, but such experiments are usually time-consuming and expensive to perform. In addition, limitations exist when extrapolating data obtained from a scaled model to a prototype, such as scaling, friction, and viscosity; experimental errors are also inevitable. Accurate definition of initial conditions, mapping of coastal features. Ocean data collected by remote sensing can be used to better understand the oceans, to provide knowledge to best manage ocean resources, and to provide useful, accurate information for commercial and recreational applications. CODAR is a land-based HF remotely sensing radar system which measures the near-surface ocean currents in a coastal area with fine temporal and spatial resolutions. One such system, consisting of two radar masts, was deployed on Galway Bay (see Figure 1) in 2011. This system is capable of monitoring the surface currents and wave parameters over most of inner Galway Bay. Measurements obtained from the CODAR system are in near real time. When a radar signal scatters off a wave whose wave length is exactly equal to half of the transmitted signal wavelength, the radar signal can return measurement information to the radar receiver [21][22][23]. A single HF radar station determines radial components of surface currents toward and away from that station. Total surface currents are computed and displayed as vector fields, by combining the radial surface current velocity components from two or more different masts [24,25].
Shulman and Paduan [26] assimilated 33-h low-pass-filtered radar data and unfiltered radar data into a model for the Monterey Bay area using OI; they found good agreement with moored current observations. Moreover, they assimilated radial and total surface currents separately. Results indicated that assimilation of radial surface currents extended the range of influence of the data into regions covered by only one HF radar site. Since the authors focused on enhancing the surface flow fields in the whole inner Galway Bay area, the unfiltered total vector fields with larger spatial coverage were employed to build a data assimilation forecasting system using nudging data assimilation technique in this work.
The CODAR system provides rich datasets in time and space. They can be used to explore the dynamical process of surface currents, to implement data assimilation and to validate numerical models [27][28][29]. The temporal and spatial resolutions of surface currents in the Galway Bay domain are sixty minutes and 300 m, respectively. The operating frequency of both radar stations in Galway Bay at C1 and C2, as shown in Figure 1, is 25 MHz. Advantages of remote sensing CODAR system over the conventional techniques include: synoptic coverage, repeated observations and high temporal and spatial resolution. Shulman and Paduan [26] assimilated 33-h low-pass-filtered radar data and unfiltered radar data into a model for the Monterey Bay area using OI; they found good agreement with moored current observations. Moreover, they assimilated radial and total surface currents separately. Results indicated that assimilation of radial surface currents extended the range of influence of the data into regions covered by only one HF radar site. Since the authors focused on enhancing the surface flow fields in the whole inner Galway Bay area, the unfiltered total vector fields with larger spatial coverage were employed to build a data assimilation forecasting system using nudging data assimilation technique in this work.
The CODAR system provides rich datasets in time and space. They can be used to explore the dynamical process of surface currents, to implement data assimilation and to validate numerical models [27][28][29]. The temporal and spatial resolutions of surface currents in the Galway Bay domain are sixty minutes and 300 m, respectively. The operating frequency of both radar stations in Galway Bay at C1 and C2, as shown in Figure 1, is 25 MHz. Advantages of remote sensing CODAR system over the conventional techniques include: synoptic coverage, repeated observations and high temporal and spatial resolution.
Although data transfer delay usually exists when processing radar results, it is less than one hour in the case of the Galway Bay CODAR system. In addition, radar data were not assimilated into models during forecasting periods. Once the data assimilation model was developed, the model can be employed to produce actual forecasts.
Measurements of surface currents in Galway Bay have been validated with independent observations from Acoustic Doppler Current Profile (ADCP) by O'Donncha et al. [20], and Ren et al. [13].

Numerical Model
The numerical model, EFDC, was used to simulate the hydrodynamic circulation of Galway Bay. EFDC was developed at the Virginia Institute of Marine Science by the U.S. Environmental Protection Agency (EPA). It comprises four linked modules: hydrodynamic, 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 motion 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 timestepping [30][31][32]. The model uses a sigma vertical coordinate system, and either regular, or curvilinear orthogonal horizontal coordinates. The model had been applied to a variety of modelling studies of rivers, lakes, estuaries, and coastal regions [33][34][35].
The dynamics within Galway Bay are dominated by oceanic flows into the bay from the Atlantic shelf and wind driven currents. Oceanic flows enter and exit the bay mainly through the sounds around the Aran Islands. Meteorological parameters are strongly affected by the Atlantic weather system [36]. The main wind direction in this area is southwest [37]. In this research, a barotropic model of Galway Bay (see Figure 1) was developed using a regular grid coordinate system; a 150 m horizontal spatial resolution was employed, yielding 380 × 241 grid cells. 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. A detailed description on setting up vertical layer structure for the Galway Bay can be found in the study by Ren et al. [38]. The meteorological forcing data (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. 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). Tidal water elevation time series generated from Oregon State University Tidal Inversion Software (OTIS) were used to define the tidal forcing on the western and southern open boundaries (see Figure 1) in the model [39,40].
The simulation period was from Julian Day 211 to 230 in 2013; this period was subdivided as follows: (i) spin-up period, Julian Day 211-220; (ii) data assimilation period, Julian Day 220-228 01:00; (iii) forecasting period, after Julian Day 228 01:00. During this simulation period, the measured radar data have high density in time and space.

Data Assimilation
Nudging is a sequential data assimilation algorithm that combines model background states with measurement states in a linear formula. Model background states were calculated from surface velocity components from the EFDC model. A nudging term was introduced into the equations of Remote Sens. 2017, 9, 932 6 of 22 motion using the difference between model background states and observation states [5]. The analysis equation has been conceptually expressed as [41]: where U is the model background states; U 0 is the observation states from radars; λ is the nudging parameter; (physics) denotes the physical process mathematically described in the numerical model. The nudging parameter is defined in the following empirical equation [41]: where r is the distance between model grid point and the observation data location; (t − t 0 ) is the difference between assimilation and observation time; t a is an assimilation timescale, which determines the strength of the nudging parameter; t d is a damping time scale for the nudging term; R n is a nudging is the exponential decay parameter, which controls the depth of influence of the nudging parameter; z is the water depth (m); z d is the depth of influence (m). Nudging algorithms are approximations of the OI data assimilation scheme; the main differences lie in that the Kalman gain used in the OI algorithm is substituted by the nudging parameter λ.
In this research, surface currents measured by the HF radar system were assimilated into the EFDC model using the nudging data assimilation algorithm. Houser et al. [42] found that better results were obtained through nudging towards a gridded analysis using spatially interpolated precipitation data into a hydrologic model, compared to others. Rusu and Guedes Soares [43] found that the accuracy of predictions was enhanced as the amount of the data assimilated increased. Moreover, in order to ensure consistency between model background states and measurements in space during the data assimilation process, nudging towards a gridded analysis presented by Stauffer and Seaman [44] was applied in this research through spatially bilinearly interpolating radar data onto model grid points. Because the spatially interpolated HF radar measurements were collocated with the model grid point and were linearly interpolated over time to every assimilation timestep, e measurements. Thus, the nudging parameter λ can be simplified to: The nudging parameter λ was set to zero at all grid points where HF radar observations were not available. The above expression of nudging parameter, λ, was used by Gopalakrishnan and Blumberg [5].

Implementation of Data Assimilation
A complete data assimilation system is comprised of three parts: a numerical model, a data assimilation algorithm, and observational data sets. Since the authors focused on simulating the surface currents in Galway Bay, in this research, a subroutine including a nudging data assimilation algorithm was encoded and combined into the EFDC model.
Since surface current vectors are the parameters of primary interest, an overview of the procedure for calculating surface currents in EFDC is shown in Figure 2. A complete data assimilation system is comprised of three parts: a numerical model, a data assimilation algorithm, and observational data sets. Since the authors focused on simulating the surface currents in Galway Bay, in this research, a subroutine including a nudging data assimilation algorithm was encoded and combined into the EFDC model.
Since surface current vectors are the parameters of primary interest, an overview of the procedure for calculating surface currents in EFDC is shown in Figure 2.  updating of surface flow fields via assimilating radar data using nudging data assimilation was carried out during the process of solving the hydrodynamic equations. HDMT.FOR is the main routine in EFDC used to calculate surface currents, so the assimilation process of radar measurements was incorporated into this routine. However, the hydrodynamic computation within EFDC is quite complex. In order to better show the updating of the surface velocity vectors via assimilating the radar data in the model, the hydrodynamic and mass transport calculation processes in the HDMT.FOR routine was extended and details are presented in Figure 3. Updating of surface flow fields via assimilating radar data using nudging data assimilation was carried out during the process of solving the hydrodynamic equations. HDMT.FOR is the main routine in EFDC used to calculate surface currents, so the assimilation process of radar measurements was incorporated into this routine. However, the hydrodynamic computation within EFDC is quite complex. In order to better show the updating of the surface velocity vectors via assimilating the radar data in the model, the hydrodynamic and mass transport calculation processes in the HDMT.FOR routine was extended and details are presented in Figure 3.  To assimilate the measured surface current data into the model, the shaded subroutine, as shown in Figure 3, was added. Implementation of the data assimilation procedure is presented in Figure 4.
Interpolation of radar data (step 1) and calculation of weight factor (λ) (step 2) were conducted before undertaking data assimilation. Updating of model background states was carried out at each assimilation timestep.  To assimilate the measured surface current data into the model, the shaded subroutine, as shown in Figure 3, was added. Implementation of the data assimilation procedure is presented in Figure 4. Interpolation of radar data (step 1) and calculation of weight factor (λ) (step 2) were conducted before undertaking data assimilation. Updating of model background states was carried out at each assimilation timestep.

Sensitivity Experiments
In order to develop a good data assimilation forecasting model for the Galway Bay domain using the HF radar data, the procedure can be described as a twofold process: firstly, appropriate nudging parameters (assimilation time scale t and the depth of influence z ) were examined, and appropriate values defined based on RMSE values between surface currents and HF radar data during the hindcasting period. Secondly, sensitivity tests of model forecasts of data assimilation cycle lengths' influence were undertaken. The authors are not aware of any other such sensitivity analyses that had been carried out to improve model forecasts using nudging data assimilation technique.

Tests of Nudging Parameters
In order to produce satisfactory results when measurements were assimilated into models using the nudging algorithm, different researchers have used different values of the nudging parameters. For example, Fan et al. [16] set the assimilation scale t = 675 s and the depth of influence of data assimilation as z = 10 m when assimilating drift and satellite data. Lin et al. [18] used the same values of nudging parameters as Fan et al. [16] in their data assimilation system; Gopalakrishnan and Blumberg [5] set values at 1800 s and 2 m, respectively. In this research, assimilation scale and depth of influence were optimized based on more than twenty representative numerical experiments, as presented in Table 1. This was the first time sensitivity analysis had been performed on the parameters of a nudging data assimilation algorithm using radar data. Twenty versions of data assimilation models were run using different values for the nudging parameters. The radar surface currents were assimilated hourly into these test models. Model NDA0, the "free run" with no data assimilated into this model, was used as a benchmark for comparisons.
RMSE between the radar data and the model results were calculated using Equations (4) and (6). RMSE of surface velocity east-west components (u) was firstly calculated with Equations (4) and (5), which average the values of RMSE in space and time; the same equations were used for north-south surface velocity components (v); and the total RMSE (u, v) was finally computed using Equation (6). RMSE (u, v) was used to assess the degree of agreement between modeled results and the HF radar data.

Sensitivity Experiments
In order to develop a good data assimilation forecasting model for the Galway Bay domain using the HF radar data, the procedure can be described as a twofold process: firstly, appropriate nudging parameters (assimilation time scale t a and the depth of influence z d ) were examined, and appropriate values defined based on RMSE values between surface currents and HF radar data during the hindcasting period. Secondly, sensitivity tests of model forecasts of data assimilation cycle lengths' influence were undertaken. The authors are not aware of any other such sensitivity analyses that had been carried out to improve model forecasts using nudging data assimilation technique.

Tests of Nudging Parameters
In order to produce satisfactory results when measurements were assimilated into models using the nudging algorithm, different researchers have used different values of the nudging parameters. For example, Fan et al. [16] set the assimilation scale t a = 675 s and the depth of influence of data assimilation as z d = 10 m when assimilating drift and satellite data. Lin et al. [18] used the same values of nudging parameters as Fan et al. [16] in their data assimilation system; Gopalakrishnan and Blumberg [5] set values at 1800 s and 2 m, respectively. In this research, assimilation scale and depth of influence were optimized based on more than twenty representative numerical experiments, as presented in Table 1. This was the first time sensitivity analysis had been performed on the parameters of a nudging data assimilation algorithm using radar data. Twenty versions of data assimilation models were run using different values for the nudging parameters. The radar surface currents were assimilated hourly into these test models. Model NDA0, the "free run" with no data assimilated into this model, was used as a benchmark for comparisons. RMSE between the radar data and the model results were calculated using Equations (4) and (6). RMSE of surface velocity east-west components (u) was firstly calculated with Equations (4) and (5), which average the values of RMSE in space and time; the same equations were used for north-south surface velocity components (v); and the total RMSE (u, v) was finally computed using Equation (6). RMSE (u, v) was used to assess the degree of agreement between modeled results and the HF radar data.
where u CODAR is the CODAR HF radar measured surface velocity east-west component (cm/s); u model is the modelled surface velocity east-west component (cm/s); RMSE j (u) is the RMSE of east-west velocity component in space covered by the high frequency radar system at timestep j (cm/s); N c is the number of times using Equation (4); N j is the number of calculation points at timestep j; RMSE(u) is the averaged value of east-west velocity component in space during the data assimilation period (cm/s). Values of RMSE(u), RMSE(v) and RMSE(u, v) over the data assimilation period are presented in Table 1.
Typical variations of the RMSE values in models are presented in Table 1. Other models were examined by the authors as well, but the difference of RMSE values was very small when z d was greater than 4 m. Averaged RMSE(u, v) between the nudging data assimilation test models, and the HF radar measurements over the data assimilation period, as shown in Table 1, indicated that model results were more sensitive to depth of influence z d , than the assimilation time scale t a . Table 1 shows that model NDA10, having the minimum RMSE(u, v) of 10.39 cm/s, generated the best forecast results.
It was assumed that good agreements between measurements and modeled results during assimilation period had positive impacts on forecasting. Thus, data assimilation scale t a = 1800 s and depth of influence of data assimilation as z d = 4 m were taken as the best nudging parameters.
To examine the performance of models using the nudging data assimilation algorithm, surface flow fields over the data assimilation period from four representative assimilation models presented in Table 1 were compared with model NDA0 and the radar data in Figure 5.
Remote Sens. 2017, 9, 932 10 of 22  Figure 5 shows that results from a representative model NDA9 (see Figure 5b) significantly deviated from the HF radar observations; similar deviation occurred when z is less than 3m. This indicated that inappropriate nudging parameters lead to model deviation or even instability. Patterns of surface flow fields in model NDA0, NDA11, and NDA12 were quite similar. Based on averaged surface flow fields during assimilation period, model NDA10 (Figure 5c) can yield closer results to radar data (Figure 5f). In the following studies, the values t = 1800 s and z = 4 m were used in nudging data assimilation models.  Figure 5 shows that results from a representative model NDA9 (see Figure 5b) significantly deviated from the HF radar observations; similar deviation occurred when z d is less than 3 m. This indicated that inappropriate nudging parameters lead to model deviation or even instability. Patterns of surface flow fields in model NDA0, NDA11, and NDA12 were quite similar. Based on averaged surface flow fields during assimilation period, model NDA10 (Figure 5c) can yield closer results to radar data (Figure 5f). In the following studies, the values t a = 1800 s and z d = 4 m were used in nudging data assimilation models.

Tests of Data Assimilation Cycle Lengths
The temporal resolution of HF radar measurements for Galway Bay domain is sixty minutes.
To study the influences of variations in data assimilation cycle lengths on forecasts, and to implement shorter updating of model background states, HF radar data were temporally linearly interpolated onto intervals less than sixty minutes. The reasons for using temporally interpolated radar data are, firstly, Sun [45] suggested that it was desirable to have rapid update cycles to capture the temporal change and to obtain useful information from the model background states. The radar data were not regarded as a single snap-shot by Sun [45], but each point in the three-dimensional data volume was assimilated at the timestep closest to the measurement time in a sequential scan. Secondly, in operational data assimilation forecasting systems, measurements from different resources at different timesteps can be assimilated into models. Thus, assimilation intervals in these models do not have to be at regular time intervals. Application of the temporally interpolated radar measurements is similar to a real forecasting data assimilation case using measurements from different sources. Thirdly, frequent updating using temporally interpolated radar data examines whether or not the three-dimensional coastal EFDC model is sufficiently robust for developing an operational data assimilation model forecasting system. Finally, since hourly output of surface current vectors were obtained by averaging/merging data over a specified time period in the CODAR radar monitoring system, temporal linear interpolation of the output data was somewhat analogous to an inverse process, which conveyed the averaged information to successive timesteps over the measurement period [46]. Linearly interpolated surface velocity data over short periods were assumed to more accurately describe the dynamics of the surface currents.
Since the ultimate goal of using data assimilation techniques was to enhance model performance through optimizing the integration of the measured data, the idea of using linearly interpolated measurements in time, to frequently update the model background states, was a numerical experiment designed to assess the accuracy of forecasted states. Additionally, Ren et al. [47] used constant pseudo measurements to test the influence of data assimilation cycle lengths in the same domain. The results showed that improved surface current forecasting resulted from the assimilation model using shorter updating intervals. It was the first time that linearly interpolated coastal HF surface currents were assimilated into a model to examine modelling performance. Temporal interpolation of radar data was implemented at each observational point where two successive hourly radar measurements were available. Here, five versions of the nudging data assimilation model were developed to assimilate temporally interpolated radar data with different data assimilation cycle lengths: (a) each model computational timestep (MS), (b) one minute, (c) five minutes, (d) fifteen minutes, and (e) sixty minutes. RMSE values were computed and averaged during the +12 h forecasting period (02:00-13:00 Julian Day 228) using Equations (4) and (6). Detailed description of these nudging data assimilation models is given in Table 2. All of the cycle length test models had the same initial and boundary conditions as model NDA10.  Table 2 shows that model NDA21, which updated the model background states at each model computational timestep using the temporally interpolated radar data, showed significant improvements during forecasting period compared with results from model NDA0. Meanwhile, other models updating model background states with longer data assimilation cycle lengths did not change much with respect to forecasting performance in models. Values of RMSE for both surface velocity components in model NDA21 were less than other models with longer data assimilation cycle length. Improvement of averaged RMSE was 10% and 15% for east-west and north-south velocity components separately during a 12-h forecasting period, in comparison with model NDA0.

Forecast Assessments of Assimilation Cycle Lengths
The above tests showed that assimilating HF radar data at each model computational timestep can improve surface flow fields using a nudging algorithm. In order to quantitatively and qualitatively compare and assess assimilation model results with the "free run" and radar data, mean surface flow fields from assimilation cycle length test models, "free run" and radar data are presented in Section 4.1. A statistical assessment method, namely DASS, of each surface flow field during a six-hour forecasting period based on MSE, is computed and presented in Section 4.2. Section 4.3 presents an analysis of the correlation of AKE between model results and radar data over a two-day forecasting period. Time series of both surface velocity components' forecasts from models at one assimilation point are also compared with radar data in Section 4.4.

Assessment of Mean Surface Flow Fields
An implicit assumption when using data assimilation techniques in numerical models was that a combination of measured data with model background states will improve modelling performance, especially during forecasting periods. In this work, we were interested in enhancing the forecasting of surface flow fields in Galway Bay through assimilating the radar data. In order to explore improvements using a nudging data assimilation algorithm over the domain during a forecast period, averaged surface flow fields during a +12 h forecasting period (start from Julian Day 228, 2013 02:00) were calculated. All nudging data assimilation models had the same conditions, except data assimilation cycle lengths as presented in Table 2. The results are presented in Figure 6, comparing different forecasts against radar data.
Averaged surface flow fields over the half-day forecasting period Julian Day 228 02:00-13:00 show visually that model NDA21 (Figure 6b), updating at each model computational timestep, produced closer surface flow fields to the HF radar data (see Figure 6g). Other nudging models (NDA10, NDA22-NDA24) using longer data assimilation cycle lengths did not show significant improvements compared with model NDA0 and the radar data. RMSE values of averaged velocity components over the half-day forecasting between the radar data and the nudging data assimilation models were calculated. RMSE from model NDA21 was smaller than model NDA0. RMSE values were 7.36 cm/s and 5.21 cm/s for the east-west and north-south velocity component of model NDA21, respectively. Improvement of averaged total RMSE (u, v) was 12% during the half-day forecasting, in comparison with the model NDA0. However, there were no distinct improvements from the other nudging data assimilation models using longer assimilation cycle lengths.
In summary, averaged surface flow fields and values of RMSE showed that model NDA21 used to update model background states at each model computational timestep significantly improved the model forecasting. Thus, temporal linear interpolation was significant and useful in this context. NDA21, respectively. Improvement of averaged total RMSE (u, v) was 12% during the half-day forecasting, in comparison with the model NDA0. However, there were no distinct improvements from the other nudging data assimilation models using longer assimilation cycle lengths.
In summary, averaged surface flow fields and values of RMSE showed that model NDA21 used to update model background states at each model computational timestep significantly improved the model forecasting. Thus, temporal linear interpolation was significant and useful in this context.

Data Assimilation Skill Score Assessment
In order to quantitatively and concisely assess the improvements of using nudging data assimilation, a DASS based on MSE was firstly calculated over time at every HF radar measurement point, and then averaged across the data assimilation domain to evaluate modelling performance.

Data Assimilation Skill Score Assessment
In order to quantitatively and concisely assess the improvements of using nudging data assimilation, a DASS based on MSE was firstly calculated over time at every HF radar measurement point, and then averaged across the data assimilation domain to evaluate modelling performance. The DASS for the north-south component of surface currents can be expressed as [48,49]: where v DA and v 0 are the north-south velocity component from the model with and without data assimilation, respectively; v CODAR is the measured north-south velocity component from the HF radar system. An analogous version of Equation (7) was also used to calculate DASS EW for surface east-west velocity component. If DASS is greater than zero, then it indicates that the data assimilation model improved the forecasting compared with model NDA0. If DASS is less than zero, then it means the data assimilation contaminated the basic dynamic processes in the model resulting in deterioration in model accuracy. The time series of DASS for both surface current components are shown in Figure 7.   Figure 7 shows that the DASS values for both surface current components were greater than zero over a six-hour forecasting period. This indicated that the application of nudging data assimilation to combine radar data into EFDC model had positive influences on forecasting. During the first three hours (02:00-04:00), DASS values were larger than DASS . This indicated that data assimilation improved the performance of north-south surface velocity component better than east-west surface velocity component. However, during the final three hours (05:00-07:00), improvement of east-west surface velocity component exceed north-south surface velocity component. In general, the range of DASS was larger than that of DASS ; magnitudes of DASS decreased over time, while magnitudes of DASS increased over the first six hours, and then decreased. This was due to that model adaptability in computing east-west surface velocity component after implementing data assimilation procedure was different from that of north-south surface velocity component. A similar trend of skill score to the DASS was obtained when assessing performance of a sea-ice model after implementing data assimilation by Levy, et al. [50]. Moreover, a similar trend of weighted skill score existed in a global weather data assimilation system by Clayton et al. [9]. Averaged DASS in the data assimilation domain during the +6 h forecast was 26%; the averaged DASS was 33%. The analysis indicates that both surface velocity components were significantly improved during the +6 h forecasting when using the nudging data assimilation algorithm to update the model background states at each model computational timestep.
So, although counterintuitive, it is not unusual for DASS to increase with time during certain forecasts. In this case, it is likely that the increase in DASS is due to the relative changes in wind and tidal forcing functions.

Averaaged Kinetic Energy Assessment
Assessment of model results based on AKE values is a convenient way to determine fundamental correlation between assimilation model results and observations. Breivik and Satra [15], and Xu et al. [51], previously applied AKE to assess model results. The spatially averaged AKE  Figure 7 shows that the DASS values for both surface current components were greater than zero over a six-hour forecasting period. This indicated that the application of nudging data assimilation to combine radar data into EFDC model had positive influences on forecasting. During the first three hours (02:00-04:00), DASS NS values were larger than DASS EW . This indicated that data assimilation improved the performance of north-south surface velocity component better than east-west surface velocity component. However, during the final three hours (05:00-07:00), improvement of east-west surface velocity component exceed north-south surface velocity component. In general, the range of DASS NS was larger than that of DASS EW ; magnitudes of DASS NS decreased over time, while magnitudes of DASS EW increased over the first six hours, and then decreased. This was due to that model adaptability in computing east-west surface velocity component after implementing data assimilation procedure was different from that of north-south surface velocity component. A similar trend of skill score to the DASS EW was obtained when assessing performance of a sea-ice model after implementing data assimilation by Levy, et al. [50]. Moreover, a similar trend of weighted skill score existed in a global weather data assimilation system by Clayton et al. [9]. Averaged DASS EW in the data assimilation domain during the +6 h forecast was 26%; the averaged DASS NS was 33%. The analysis indicates that both surface velocity components were significantly improved during the +6 h forecasting when using the nudging data assimilation algorithm to update the model background states at each model computational timestep.
So, although counterintuitive, it is not unusual for DASS to increase with time during certain forecasts. In this case, it is likely that the increase in DASS is due to the relative changes in wind and tidal forcing functions.

Averaaged Kinetic Energy Assessment
Assessment of model results based on AKE values is a convenient way to determine fundamental correlation between assimilation model results and observations. Breivik and Satra [15], and Xu et al. [51], previously applied AKE to assess model results. The spatially averaged AKE (E k ) across the study domain was calculated to assess improvements when using nudging data assimilation. E k was calculated as follows [15,49,52]: where A is the domain covered by the CODAR system; ds is the integral cell area; U is total surface speed. Water density is assumed constant, and therefore, it is omitted from Equation (8).
The AKE from models and observations was computed in the data assimilation domain at each model forecasting timestep. Hourly data comprised of 48 points from each dataset over a 2-day forecasting period (start from Julian Day 228 02:00) were used for this analysis, which was sufficient to ensure the analysis is meaningful. Time series of AKE values were calculated using Equation (8).
Correlation of AKE between model NDA21 and radar data was improved by 12% in comparison with the model NDA0.

Assessment of Surface Velocity Components
To further compare the improvements among models and radar data in time, time series of surface velocity components at location A (see Figure 1) from six representative nudging data assimilation models are respectively shown in Figures 8 and 9. where A is the domain covered by the CODAR system; ds is the integral cell area; U is total surface speed.
Water density is assumed constant, and therefore, it is omitted from Equation (8).
The AKE from models and observations was computed in the data assimilation domain at each model forecasting timestep. Hourly data comprised of 48 points from each dataset over a 2-day forecasting period (start from Julian Day 228 02:00) were used for this analysis, which was sufficient to ensure the analysis is meaningful. Time series of AKE values were calculated using Equation (8).
Correlation of AKE between model NDA21 and radar data was improved by 12% in comparison with the model NDA0.

Assessment of Surface Velocity Components
To further compare the improvements among models and radar data in time, time series of surface velocity components at location A (see Figure 1) from six representative nudging data assimilation models are respectively shown in Figures 8 and 9. Figures 8 and 9 show significant improvements over time of both velocity components appeared in the best data assimilation model NDA21 compared with model NDA0 at point A. Other nudging data assimilation models (NDA10, NDA22-NDA24), with longer data assimilation cycle lengths, had quite a similar trend to model NDA0 during the forecasting period. Although the trend of the north-south velocity component from model NDA21 was different from the radar data (shown in Figure 9), its trend was closer to the radar data than results from the other models. Moreover, the magnitudes of the north-south velocity component were quite small during the analysis period. Because the north-south velocity component was mainly affected by wind forces, using constant wind data in time and space in the model may lead to errors.   Figures 8 and 9 show significant improvements over time of both velocity components appeared in the best data assimilation model NDA21 compared with model NDA0 at point A. Other nudging data assimilation models (NDA10, NDA22-NDA24), with longer data assimilation cycle lengths, had quite a similar trend to model NDA0 during the forecasting period. Although the trend of the north-south velocity component from model NDA21 was different from the radar data (shown in Figure 9), its trend was closer to the radar data than results from the other models. Moreover, the magnitudes of the north-south velocity component were quite small during the analysis period. Because the north-south velocity component was mainly affected by wind forces, using constant wind data in time and space in the model may lead to errors.

Assessment of Forecasted Surface Flow Fields
Although the above mean surface flow fields showed surface flow patterns during a half-day period, surface flow fields at a representative forecasting time, 03:00 Julian Day 228, are presented in Figure 10, to further illustrate model performance after implementing data assimilation.

Assessment of Forecasted Surface Flow Fields
Although the above mean surface flow fields showed surface flow patterns during a half-day period, surface flow fields at a representative forecasting time, 03:00 Julian Day 228, are presented in Figure 10, to further illustrate model performance after implementing data assimilation.

Assessment of Forecasted Surface Flow Fields
Although the above mean surface flow fields showed surface flow patterns during a half-day period, surface flow fields at a representative forecasting time, 03:00 Julian Day 228, are presented in Figure 10, to further illustrate model performance after implementing data assimilation.  The east-west trend of surface currents was stronger in model NDA21 than model NDA0. The general trend in model NDA21 was closer to the radar measured flow field. Improvements of vector direction of surface currents over the data assimilation domain in model NDA21 were 20% compared with model NDA0 at 3:00. This further demonstrated that updating the model background states at each model computational timestep, using nudging assimilation algorithm, significantly enhanced forecasting of surface flows.

Discussion
As the difference between observations and model results was weighted using the factor λ in the nudging data assimilation technique, and then used to update model background states, the complex definition of observation error covariance and model error covariance used in OI and 4D-VAR techniques was not required. Thus, the application of the nudging algorithm is simpler and more efficient than others, such as OI and 4D-VAR data assimilation schemes. This is of great importance for a real-time operational data assimilation forecasting model, which requires timely prediction information for practical operations, such as search and rescue, and oil spill response. Selection of the most appropriate data assimilation algorithm is important when exploring the potential for using HF radar data to improve model performance.
According to a DASS forecast skill score, improvement of velocity component prediction from a nudging data assimilation model updating at each model computational timestep, was more significant than the assimilation model by Ren et al. [13], which indirectly corrected wind stress using HF radar data for the same research domain. Averaged DASS values over six-hour forecasts obtained by Ren et al. [13] were 1% and 23% for east-west and north-south velocity components, respectively; DASS of 26% and 33% were obtained for velocity east-west and north-south components, respectively, in this work, during the same prediction period. Although indirect data assimilation via correcting wind stress using HF radar data significantly improved forecasting of velocity of the north-south component, the influence of data assimilation on velocity of the east-west component was not as good as using a nudging algorithm. This indicated that the approach proposed in this work was more suitable for producing accurate patterns of surface flow fields for the Galway Bay area. Additionally, Gopalakrishnan and Blumberg [5] used a similar nudging data assimilation scheme to combine HF radar data with model background states; improvements of forecasted velocity east-west (or north-south) components in this work were more significant than their assimilation model results based on DASS values. DASS in Galway Bay for velocity east-west (or north-south) component was 26% (33%), in comparison with their values at 18% (7%). This was because the data assimilation cycle length was small, i.e., at each model computational timestep in this work. In short, frequent updating of model background states using a nudging algorithm positively constrains circulation of surface flow fields during the forecasting period, than those assimilation models with longer data assimilation cycle lengths.
The best data assimilation model (NDA21) had a much smaller RMSE value than RMSE values (range 10~20 cm/s) obtained by Zhao et al. [53] using quasi-ensemble Kalman filter with Canadian quick covariance. This may be because hourly radar data were assimilated into models in their work, but a shorter data assimilation cycle length was applied herein. It also showed that application of quasi-ensemble Kalman filter for assimilating radar data was more challenging than using a nudging data assimilation algorithm.
According to the AKE assessment criteria, improvement of AKE over a 2-day forecasting period of 12% in the best nudging data assimilation model in this work, was comparable to the AKE value of 16% over a 2-h forecast obtained by Breivik and Satra [15]. However, at 4-h and 6-h forecasts, correlation of AKE between radar data and assimilation model significantly decreased from 0.58 to 0.38 in Breivik and Satra [15]. Multiple assessment criteria showed that developed nudging data assimilation updating at each model computational timestep using radar data produced comparable, or better results than others. In addition, sensitivity experiments on nudging parameters presented in Table 1 indicated that hindcasting results were more sensitive to the influence depth than the assimilation time scale. Water depths of the assimilation domain covered by CODAR system ranged from~5-35 m. To examine model performance at spatial locations with different water depths, relationships between depth of influence and the value of e [ z z d ] are presented in Table 3.  Table 3 shows that the nudging weighing factor, λ, is equal to 1 when water depth is 30 m. Magnitudes of the nudging weighing factor, λ, had the same trend as water depth varies. Significant deviations in Figure 5b existed at some points, due to incorrectly assimilating the increment U 0 − U to model background states using Equation (1). Exploration of the sensitivity of influence depth over a domain needs to be carried out in future research.
As presented in Table 2, value of the nudging parameter λ, determined by z, z d , and t a , was different at each model grid. The nudging parameter λ had similar trends relative to z and z d , when t a was set different values. In order to further investigate relationships among them, an example as t a = 1800 s is shown in Figure 11. To clearly show the main variation range of λ relative to z, z d and t a , only parts of plots were shown in Figure 11. The larger the value of z, then the larger λ was, given a fixed value of z d . It indicated that those points with deeper water depth were nudging to a larger extent. This lead to significant deviation of analysis states from radar data. Use of varying z d in space can deal with overly nudging at those points, with deeper water depth z. Moreover, amount of λ, whose value was greater than 10, increased as z d decreased. It indicated that z d needed to be large enough, such as 4 m, as in this work. This can ensure that model background states at fewer points were overly nudged than those when z d was smaller. Thus, presented sensitivity experiments in this work were a must for other cases when using similar nudging data assimilation algorithms.
Remote Sens. 2017, 9, 932 18 of 22 by CODAR system ranged from ~5-35 m. To examine model performance at spatial locations with different water depths, relationships between depth of influence and the value of e are presented in Table 3.  Table 3 shows that the nudging weighing factor, λ, is equal to 1 when water depth is 30 m. Magnitudes of the nudging weighing factor, λ, had the same trend as water depth varies. Significant deviations in Figure 5b existed at some points, due to incorrectly assimilating the increment (U − U) to model background states using Equation (1). Exploration of the sensitivity of influence depth over a domain needs to be carried out in future research.
As presented in Table 2, value of the nudging parameter λ, determined by z, z , and t , was different at each model grid. The nudging parameter λ had similar trends relative to z and z , when t was set different values. In order to further investigate relationships among them, an example as t = 1800 s is shown in Figure 11. To clearly show the main variation range of λ relative to z, z and t , only parts of plots were shown in Figure 11. The larger the value of z, then the larger λ was, given a fixed value of z . It indicated that those points with deeper water depth were nudging to a larger extent. This lead to significant deviation of analysis states from radar data. Use of varying z in space can deal with overly nudging at those points, with deeper water depth z. Moreover, amount of λ, whose value was greater than 10, increased as z decreased. It indicated that z needed to be large enough, such as 4 m, as in this work. This can ensure that model background states at fewer points were overly nudged than those when z was smaller. Thus, presented sensitivity experiments in this work were a must for other cases when using similar nudging data assimilation algorithms. In addition, models using the same parameters as in model NDA9, NDA11, and NDA12, were performed to assimilate the interpolated radar data at each model computational timestep. On the basis of computation and comparison of RMSE values between assimilation model results and In addition, models using the same parameters as in model NDA9, NDA11, and NDA12, were performed to assimilate the interpolated radar data at each model computational timestep. On the basis of computation and comparison of RMSE values between assimilation model results and radar data during a 12 h forecasting period, the model using the same nudging parameters as in model NDA9 to assimilate radar data at each model computational timestep overly nudged the analysis states, leading to model instability, and forecasting was significantly larger than radar data. Models using the same nudging parameters as in model NDA11 and NDA12 were developed to assimilate the interpolated radar data at each model computational timestep, respectively. Results from these models were quite similar to outputs from the best model NDA21, but model NDA21 was slightly better than others, in terms of RMSE values between model outputs and radar data during a 12 h forecasting period. Experiments performed in this work indicated that the same nudging parameters can be used to develop nudging data assimilation models with different cycle lengths.

Conclusions
Real-time remotely sensing HF radar measurements were blended into a numerical model EFDC using a nudging data assimilation algorithm. The most appropriate nudging parameters were determined based on obtaining minimum RMSE values between modeled results and radar data during a data assimilation period with hourly assimilation. Assimilation of the remotely sensing radar data at each model computational timestep resulted in a good forecasting performance compared with a "free run" model and models with longer data assimilation cycle lengths. The main conclusions from the research are: (1) Sensitivity tests performed on the nudging data assimilation parameters suggest that an influence depth of 4 m and an assimilation timescale of 1800 s were the suitable values for developing the data assimilation system for Galway Bay based on RMSE analysis during a hindcasting period. Further analysis indicated that surface flow fields were more sensitive to the depth of influence than the assimilation time scale. This is the first study that has investigated the effects of the system parameters on results, and to provide guidance for future application of this algorithm. (2) The research showed that modelling performance was also quite sensitive to data assimilation cycle lengths. Assimilation of radar data at each model computational timestep significantly improved model forecasting in comparison with using longer data assimilation cycle lengths. The RMSE of the averaged velocity components over a half-day forecasting period between the radar data and the best nudging data assimilation model was smaller than model forecasts using longer cycle lengths. The averaged east-west and north-south velocity components from the best model (NDA21) were improved by 10% and 15%, compared with the "free run", respectively. However, there were no distinct improvements in RMSE between the other models with longer data assimilation cycle lengths. The results presented herein demonstrated better performance than all other published results in investigating this area. Additionally, results indicated that the same nudging assimilation parameters can be used for models with different cycle lengths. (3) The calculation of surface flow fields were enhanced during the forecasting period when data assimilation was employed at each model computational timestep. Surface flow fields from other assimilation models with longer cycle lengths did not significantly improve compared with model NDA0. This was the first time such a sensitivity analysis was performed. (4) The values of DASS generated during this research proved that model NDA21 improved model performance by 26% and 33% for east-west and north-south velocity components, respectively, during a +6 h forecasting period for Galway Bay; these were significant improvements. (5) Analysis of AKE further proved that updating model background states at each model computational timestep resulted in forecasting improvements. Correlation of AKE between the best nudging data assimilation model NDA21 and the radar data was improved by 12% compared with model NDA0 over a two-day forecasting. The improvement, herein, was comparable with other studies.
(6) Time series graphs at a location within Galway Bay showed that implementation of data assimilation at each model step using interpolated radar data greatly enhanced both surface velocity components during forecasting period, especially for the north-south surface velocity component. Results from other assimilation models with longer data assimilation cycle lengths did not generate distinct improvements on the model forecasts; they followed the same trend as the "free run". (7) Surface flow fields at a representative measurement timestep, 03:00 Julian Day 228, show that the best assimilation model produced significantly better surface flow circulation through the domain than the "free run" did. Improvements of vector directions in the domain were in the order of 20% at this time, compared with the "free run".
In summary, the assimilation of remotely sensed HF radar data using a nudging algorithm is a powerful tool for improving model performance. For the first time, this research presented evidence of the effects of changes in nudging parameters and data assimilation cycle lengths on hydrodynamic forecasting, when applying this algorithm. The nudging algorithm has been shown to be particularly useful when updating is applied at each model computational timestep. Significant improvements in forecasting can last more than six hours after data assimilation; these forecasts can provide useful information for a variety of applications, such as search and rescue and oil spill operations. Sensitivity analysis, as described above, should be carried out for site specific models when using remotely sensed data. In order to explore effects of assimilation with high temporal model updating and to further enhance modelling performance, non-linear interpolated temporal radar data will be conducted in future work.