Constraining Energetic Slope Currents through Assimilation of High-frequency Radar Observations

Assimilation of high-frequency (HF) radar current observations and CTD hydrography is performed with the 4D-Var analysis scheme implemented in the Regional Ocean Modeling System (ROMS). We consider both an idealized case, with a baroclinic slope current in a periodic channel, and a realistic case for the coast of Vesterålen in northern Norway. In the realistic case, the results of the data assimilation are compared with independent data from acoustic profilers and surface drifters. Best results are obtained when background error correlation scales are small (10 km or less) and when the data assimilation window is short, i.e. about 1 day. Furthermore, we find that the impact of assimilating HF radar currents is generally larger than the impact of CTD hy-drography. However, combining the HF radar currents with a few hydrographic profiles gives significantly better results, which demonstrates the importance of complementing surface observations with observations of the vertical structure of the ocean.


Introduction
Skillful ocean forecasts are of key importance for many operations at sea, especially for emergency response services such as search-and-rescue and oil spill mitigation. In particular, near-surface currents are an important input to operational drift forecast models. However, the predictability of ocean currents remains a challenge due to their turbulent nature and high spatial variability, for example, associated with eddies.
Observations of the ocean surface temperature (Wentz et al., 2000;Rayner et al., 2003) and elevation (Fu et al., 1994) from satellites have become plentiful during the last decades. New satellite observations include surface salinity and currents, but the uncertainty of these products still remains too high for use in models with high horizontal resolution, i.e. on the order of 1 km. Through efforts such as the International Argo Program (Roemmich et al., 2009), observations of the subsurface ocean are increasing in number, but still remain too few to resolve the vertical and horizontal density structure of, e.g. oceanic fronts.
Advanced data assimilation (DA) techniques developed in the field of numerical weather prediction are now used to a great extent within ocean modelling. The DA schemes range from multivariate implementations of optimal interpolation, such as Ensemble Optimal Interpolation (Oke et al., 2010) in the BLUElink forecast system (Brassington et al., 2007), to ensemble Kalman filters (EnKF) (Evensen, 2003) in the TOPAZ 4 system (Sakov et al., 2012) and variational methods (Dimet and Talagrand, 1986;Courtier et al., 1994) used in UK Met Office's FOAM (Blockley et al., 2014). Kalnay et al. (2007) and Gustafsson (2007) provide an interesting discussion on the advantages of both EnKF and the timedependent variational method, 4D-Var. The latter approach is applied in this study. In 4D-Var, the aim is to minimize a cost function, which takes both misfits between model and background state and between model and observations into account. One of the main benefits of 4D-Var lies in the fact that observations are evaluated at the correct time.
Common for these forecast systems is that few, if any, current observations are assimilated, mainly owing to a lack of observations. HF (high-frequency) radars, which can sample surface currents up to 200 km offshore, provide an excellent option for mapping surface currents in coastal areas (e.g. Paduan and Washburn, 2013;Barrick et al., 1977;Chapman et al., 1997;Gurgel et al., 1999). Real-time surface cur-Published by Copernicus Publications on behalf of the European Geosciences Union. 238 A. K. Sperrevik et al.: Assimilation of high-frequency radar currents rents from HF radars can be used for monitoring and emergency applications, but also for data assimilation to improve the ocean forecasts. Previous studies (see Breivik and Saetra, 2001;Oke et al., 2002;Paduan and Shulman, 2004;Barth et al., 2008;Zhang et al., 2010, for some examples) demonstrate the potential of HF radar current observations in ocean data assimilation systems. The aim of this study is to investigate whether assimilation of current observations from a rapidly deployable HF radar system (Kjelaas and Whelan, 2011) in a high-resolution ocean model is a feasible way to improve the regional ocean forecast during, e.g. an oil spill event.
Assimilation of HF radar currents requires an ocean forecast system which is ready to use such observations as soon as they become available. The Regional Ocean Modeling System (ROMS; http://www.myroms.org), for which a sophisticated 4D-Var assimilation system, described in detail in Moore et al. (2011a, b, c), has been developed, is implemented for a region in northern Norway. The circulation pattern in this area is dominated by two northward flowing currents, the North Atlantic Current which flows along the shelf slope, and the Norwegian Coastal Current. Typical current speeds are on the order of 0.2-0.5 m s 1 . The properties of the water masses associated with these currents are in strong contrast to one another, as the coastal current strongly depends on freshwater from the Baltic Sea and the Norwegian fjords. During winter the coastal water is also noticeably colder than the saltier and warmer water masses of Atlantic origin further offshore. The region just west of Vesterålen, where the continental shelf narrows considerably, is characterized by high eddy kinetic energy levels (Isachsen et al., 2012), which makes this a challenging area for ocean predictions. Additionally, there is a strong tidal signal in the area (Moe et al., 2002), with tidal ranges of 2-3 m. Initial tests using an idealized model were made before the tests with a realistic model configuration. These tests were motivated by the wish to investigate more in detail the possible impact of HF radar observations in constraining unstable, highly energetic and narrow currents on steep slopes, such as we find offshore of Vesterålen. These idealized experiments were carried out prior to the field experiment and provided useful information about the expected impact of assimilating real observations, as well as starting points for tuning the data assimilation system with regard to values of e.g. assimilation window length and horizontal error correlation length scales.
We use the incremental, strong-constraint (perfect model assumption) IS4DVAR driver of ROMS, in which the cost function is minimized iteratively by applying a conjugate gradient algorithm, so-called inner loops. Outer loops allow for relinearization of the full model trajectory where the intermediate estimates of the cost function are taken into account, and is a way to account for nonlinearities not represented in the tangent linear and adjoint models. ROMS IS4DVAR is well documented in a series of studies, such as Powell et al. (2008), Broquet et al. (2009), Zhang et al. (2010 and Zavala-Garay et al. (2012). As a consequence of the perfect model assumption, the analyses produced by IS4DVAR corresponds to a trajectory of the forecast model. Thus, this approach is not capable of correcting discrepancies caused by flaws in the numerical model itself. We show that the skill of the predictions produced by the ocean forecast system indeed increases after just one assimilation cycle when HF radar total currents are included. Including in situ temperature and salinity profiles to the assimilated observational data set adds an additional constraint on the circulation.
We start by describing the observational data sets in Sect. 2. In Sect. 3 we present an idealized assimilation experiment, before a realistic experiment is presented and compared with independent observations in Sect. 4. A summary and some concluding remarks are given in Sect. 5.

Field campaign and observation network
An array of three SeaSonde HF radars, manufactured by CODAR Ocean Sensors, was deployed along the coast of Vesterålen in northern Norway in the first week of March 2013. All three stations were demobilized in the beginning of June 2013. The operating frequency of the radars was 13.525 MHz, thus, measuring the Bragg backscatter from waves with wavelength of about 11 m. Paduan et al. (2006) provide an excellent analysis of sources of errors in the raw radial currents. When the radials are combined to form gridded fields of total current vectors, an additional error, geometrical dilution of precision (GDOP) (Chapman et al., 1997), is introduced. Each observation data file transmitted by the HF radar system contained estimates of the observation errors, which are a combination of signal to noise rations, velocity variances within range cells, and GDOP for two or more systems. Typically, these errors were less than 20 % of the observed values within a distance of about 40 km of the radars, gradually increasing to almost 100 % near the limits of their range. The average relative error is found to be similar to what would be expected if the errors were mainly due to GDOP. As the currents observed by HF radars include tidal currents, the observations are usually subjected to a correction of the tidal signal before assimilation to account for discrepancies in observed and modelled tides, as described in Zhang et al. (2010). However, as our time series is too short to provide a good estimate of the observed tidal signal, no corrections of this kind have been made. This would, however, also be the case during a realistic event, and our results are thus a demonstration of the potential impact of a rapidly deployable HF radar system.
The main motivation for deploying the radars in this part of Norway was the co-location with the annual cod stock assessment cruise of the Institute of Marine Research (IMR) (see Fig. 1 a total of 96 CTD profiles were collected (CTD stations are shown in Fig. 1). Cross-shelf sections allow for sampling of water masses of different characteristics. The near-surface observations vary from 2 C and a salinity of 33 inside the archipelago to 6 C and a salinity of 35 at the furthest point offshore. Seven surface drifters from MetOcean, Canada, were deployed from the research vessel (R/V) Johan Hjort as it passed the area covered by the HF radars (see Fig. 2). These were iSLDMB (iridium Self Locating Datum Marker Buoy) drifters, which have a drogue centred at 65 cm below the surface. The precision of the GPS positions are approximately 10 m, implying an error circle with a radius of 10 m. The analysis in Sect. 4.3 focuses on 3 h long trajectory segments during which the buoys travelled an average of ⇠ 2700 m. The errors are thus typically less than 1 % of the observed distance.
Prior to the cod stock assessment cruise we deployed a mooring with three separate acoustic Doppler current profilers (ADCPs). The total water depth at the site was 86 m, its location is shown in Fig. 2. Two Anderaa ADCPs of 600 kHz were mounted at a depth of about 40 m, one upward looking and one downward looking, in addition to an upward looking 1 MHz Nortek Aquadopp Profiler at about 10 m depth. The observed ocean currents at 1 m depth from the Nortek Aquadopp ADCP had a variability (standard deviation) of 17 cm s 1 . The accuracy is 2 cm s 1 for 30 min averaged samples. This accuracy is obtained as the mean error of the 30 min averages, given by the standard deviation of bootstrapped samples in each averaging period. campaign is described in more detail in Röhrs et al. (2015) and Christensen et al. (2013).

Idealized experiments
Prior to applying the assimilation system to a realistic case we investigate various options associated with the IS4DVAR scheme using an idealized setup. More specifically, the IS4DVAR scheme is tested for an idealized case of strongly nonlinear, unstable baroclinic flow along a steep slope. The impact of varying the horizontal error correlation scales, which are isotropic, the number of inner and outer loops, as well as the length of the assimilation window is assessed to provide an indication on how to configure the assimilation system for the realistic case later on. The tests are set up such that the model system is essentially without any systematic error and the evaluation of the test cases is primarily based on the standard deviation of the analysis error. The main motivation of these experiments is to test the applicability of IS4DVAR for such energetic and nonlinear flow as we find in our study area. Focus is therefore on the ability of IS4DVAR to reproduce the observed flow features; hence, we have not added any errors to the synthetic observations. We use a linear equation of state, and density changes are accommodated by imposing surface and bottom heat fluxes; hence, in contrast to the realistic experiments described later on, salinity plays no role here. There are seven prognostic variables in ROMS (u, v, T , S, ⇣,ū,v) representing horizontal velocity in 240 A. K. Sperrevik et al.: Assimilation of high-frequency radar currents easterly direction, horizontal velocity in northerly direction, potential temperature, salinity, vertically averaged velocity in easterly direction, vertically averaged velocity in northerly direction, and sea surface height, respectively.

Model grid
The model domain is configured as a channel with periodic north-south boundary conditions and solid walls along the eastern and western boundaries. The grid is Cartesian with a horizontal resolution of 2.4 km. We use the f plane approximation with a constant Coriolis parameter equivalent to 65 N. The domain size is 100 ⇥ 120 interior grid points (easterly/northerly directions), and we use 35 vertical levels.
The main topographic features are (i) a shelf along the eastern boundary, with average depth of 200 m, and width of approximately 70 km, (ii) a sharp shelf break with a hyperbolic tangent profile (maximum slope is approximately 0.1), and (iii) a deep ocean floor with average depth of approximately 1800 m towards the western boundary. In order to trigger instabilities, the depth is perturbed with random values between ±50 m.

Initialization and forcing of the model
The initial conditions are uniform with T = 10 C, S = 35, ⇣ = 0 m and zero velocities. Constant heat fluxes over a period of 150 days are applied. The aim is to produce unstable baroclinic currents that will be guided by topography. To accommodate this, we apply a net bottom heat flux into the ocean on the shelf and a net surface heat flux out of the ocean over the deep ocean (Isachsen, 2011). These fluxes are exactly balanced, such that the net heat flux into the ocean is zero. Since the deep ocean is wider than the shelf, the surface fluxes are smaller than the bottom fluxes (approximately 160 and 440 W m 2 , respectively). The momentum and net freshwater fluxes are zero.
The water heated over the shelf bottom is rapidly mixed upwards through the entire water column resulting in a sharp temperature front near the shelf break. Due to geostrophic adjustment, a northward flowing, topographically controlled slope current develops. This current is baroclinically unstable and heat exchange with the deep water region is facilitated by macro turbulence (ocean eddies), see Fig. 3.
The domain is periodic so that the water masses flowing out of the domain at the northern boundary reappear at the southern boundary. As shown in Fig. 3, the mean northward surface velocity has a maximum of about 1 m s 1 over the slope, which means that a drifting object can potentially be advected through the domain in about 3 days. The flow is also characterized by baroclinic and barotropic waves and eddies with propagation speeds that are generally different from the mean flow speeds. To investigate how upstream conditions influence the dynamics in this region, we performed an adjoint sensitivity study (e.g. Moore et al., 2009; Figure 3. The mean northward surface velocity during a 3-day simulation. The region enclosed by the black line is used for assessing the performance of the data assimilation system. 2009), focusing on the advection of water masses, hence, investigating the sensitivity of the surface velocity variance to temperature in a short section over the slope in the middle of the domain. Based on our results we conclude that an upper limit of 3 days for the assimilation experiments is adequate.

Configuration of the data assimilation system
In our experiments, data assimilation is only considered for interior grid points and no adjustment of the surface/bottom fluxes or boundary conditions are made. The 4D-Var schemes implemented in ROMS also has options for multivariate background error correlations, which is based on the assumption of geostrophic balances and the baroclinic contribution to sea surface height as described in Moore et al. (2011c). In our case, the Rossby radius is small and the flow is highly unstable and energetic, implying considerable ageostrophic forcing. In addition, the main current is close to the coast and runs along very steep bathymetry; hence, we do not make use of any such options here. Correlation between state variables are, however, implicitly accounted for by the model physics. The estimates of background errors are taken from day 120 to 150 of the spin-up period, using the standard deviation of each variable in each grid point. The vertical error correlation lengths for all variables is set to 30 m. Horizontal length scales are isotropic and the same for all variables in the entire model domain (Moore et al., 2011c).

Synthetic observations
Synthetic observations are taken from a separate model simulation that has been forced with time-varying momentum and heat fluxes. An example of the difference between the simulation used for synthetic observations and the simulation that forms the basis for the assimilation experiments is shown in  . The difference between the simulations is primarily related to small-scale features associated with the baroclinically unstable current over the slope. An estimate of realistic observation errors is assigned to the synthetic observations for use in the assimilation system, but no random or systematic errors are added to the observation values themselves.
The observation locations are shown in Fig. 5. Hydrographic observations, e.g. such as those taken from a research vessel with a CTD, are simulated by taking a single vertical profile of temperature every hour, zigzagging southwards with four sections across the slope. A total of 64 temperature profiles are processed and each individual observation is assigned a constant error of 0.05 K. Two simulated HF radars are used to provide hourly total current vectors in 61 locations. These HF radar stations are positioned a distance of y a = 118 km and y b = 166 km from the southern boundary. The idealized observations are taken at positions where the beams from the two simulated HF radars intersect.
We assume that the HF radars retrieve radial currents from an effective depth, D e = 2 m, with an azimuthal resolution of 1✓ = 11.25 . Furthermore, we assume that the maximum range of the radars is R = 80 km, and that the relative observation error R associated with the radial currents is a linear function of radial distance r. The azimuthal resolution determines the number of beam directions and, combined with the maximum range R, also the number of intersecting beams from which we can estimate total current vectors.
To obtain the synthetic HF radar observation errors we have first used standard vector algebra to determine the positions were the beams intersect and then calculated the errors in the easterly and northerly directions ( GDOP ) due to geometric dilution of precision (GDOP) (Chapman et al., 1997). The total errors are then assumed given by where (u, v) are the observed speeds in the easterly and northerly directions, respectively. The errors in a real HF radar system are more complex than those given by Eq. (1), but a comparison with error statistics from a CODAR Sea-Sonde system deployed near Fedje in western Norway indicate that the values obtained from Eq. (1) are reasonable (not shown here).
Finally, in order to cover the region with the most energetic currents near the slope, the entire simulated HF radar system is translated approximately 40 km westwards from the eastern boundary, see Fig. 5.

Results of the idealized experiments
We only consider the velocities (u, v) and the temperature T in the evaluation of the idealized experiments. Furthermore, we only consider a limited region of interest similar to the HF radar coverage area (see Fig. 3) and also restrict the evaluation to the two uppermost vertical levels of the model. We evaluate a 3-day period using sequential data assimilation when the assimilation window is shorter than 72 h.
The procedure is as follows: each model simulation results in an analysis, i.e. a solution of IS4DVAR, which is compared with the "truth" as represented by the simulation that provides the synthetic observations. Average bias and standard deviation for each of the variables (u, v, T ) are calculated based on all grid points and all output times in the region of interest. Due to the idealized design of the experiments, the (u, v, T ) variables have little bias compared to what we would expect from a realistic model. The standard deviation is used when deciding which options yield the best results.  A total of 15 different model simulations with IS4DVAR have been made. In the simulations we have considered (i) the horizontal background error correlation scales, (ii) relinearization of tangent linear and adjoint models using outer loops, (iii) length of assimilation window, and (iv) the relative impact of HF radar observations compared to hydrographic observations.
The results show that using a small horizontal error correlation scale of 5 km gives best results (see Table 1); a likely explanation is that observation values are erroneously distributed across sharp gradients when larger correlation scales are used. Furthermore, using outer loops is beneficial, with no significant improvement when using more than two relinearizations. The length of the assimilation window also plays a role: slight improvement is obtained when reducing the window length from 72 to 24 h, but there is essentially no difference when the window length is further reduced to 6 h. The impact of assimilating HF radar currents is generally larger than the impact of assimilating hydrography, but in both cases we obtain improvement for the temperature. Assimilating only hydrography profiles does not improve the currents, while assimilating HF radar currents has a positive impact on both the velocities and the temperature (see Table 2).

Realistic experiments
The realistic experiments were carried out with a model application covering the region of the Lofoten and Vesterålen archipelago with 2.4 km horizontal resolution (see Fig. 6). The bathymetry of the domain resembles that of the idealized application, with a shelf in the eastern part of the domain and a steep shelf break.

Model grid and forcing data
The model domain of the realistic model simulations is a subset of the operational high-resolution ocean model at MET Norway (Albretsen et al., 2011) centred around the Lofoten and Vesterålen archipelago. The model has 35 vertical layers with increased resolution near the surface. This model setup has also been used for transport estimates of northeastern Arctic cod eggs and larvae, and a more thorough validation using in situ and drifter data is presented in Röhrs et al. (2014). As ocean data assimilation requires extensive supercomputing resources, the horizontal resolution has been decreased from 800 m to 2.4 km compared to the operational setup.
The lateral boundary conditions are retrieved from the operational setup at 800 m resolution, using 3-hourly fields of sea surface elevation, temperature, salinity, and currents. To remove fine-scale features from the high-resolution fields, that are unresolved in the coarser grid, the fields are averaged over 3⇥3 grid points before they are interpolated to the coarser grid. The simulations use open boundary conditions for sea surface elevation and barotropic currents (Chapman, 1985;Flather, 1976). For tracers and baroclinic velocities, boundary conditions as described in Marchesiello et al. (2001) are used. The adjoint and tangent linear models do not have the same options for boundary conditions, however, so that clamped boundary conditions with a sponge layer are applied for these.
As atmospheric forcing we use hourly fields of air temperature and humidity 2 m above ground, 10 m winds, sea level pressure, cloud cover, and precipitation from MET Norway's operational weather forecast at 4 km resolution (Kristiansen et al., 2009).
A spin-up model simulation was initialized at 15 February 2013 from the smoothed high-resolution fields, from the operational setup, and run through 15 April 2013. The initial conditions for the assimilation simulations were obtained from the spin-up simulation.

Configuration of the data assimilation system
The configuration of the assimilation system in the realistic experiment is based on the results from the idealized case (see Sect. 3.3). Overall, the parameters values we find to be optimal in the idealized case are also optimal in the realistic case, with the exception of the horizontal error correlation length scale, which is increased. Based on the outcome of a series of tests, we proceed with horizontal background error correlation scales of 10 km, 10 inner and 2 outer loops.
To further assess the impact of the assimilation window length on the performance of the assimilation system, a twin experiment was conducted with the realistic model setup. A model simulation covering the winter months of 2011 was used as the true ocean state. Observations of u and v, in an area similar to that of HF radar coverage during the field campaign, were collected from the true state at every hour for a duration of 72 h. Random errors, with mean value of zero, were added to the observations. A period in 2013 with similar weather conditions was identified, and used as initial conditions for the assimilation experiments. The assimilation system was run for this 72 h period several times, successively decreasing the window length from 72 to 3 h. Comparing surface currents in the resulting analyses with the "true" state, we find the assimilation window length to have little impact on analysis quality. Based on the results from the idealistic experiments, as well as results from the realistic model with real observations, a window length of 24 h has been used in the following experiments.
In our tests, we also experimented with correlated observation errors in u and v (which is inevitable due to the way total vectors are calculated). The difference between runs with correlated vs. uncorrelated errors added to the synthetic observations was marginal.
We evaluate the performance of the data assimilation system using independent observations from the field campaign (Sect. 2); that is, observations that have not been used during the data assimilation and, therefore, serve as an independent measure of skill.

Comparison with surface drifters
Ocean forecasts are used for input to trajectory models, which predict the drift of, e.g. oil spills or objects in the sea. Thus, the ability of the forecast to reproduce the trajectories of the surface drifters is a good measure of forecast skill.
To evaluate the error of the assimilated HF radar observations we compare radial HF radar currents with the corresponding component of drifter speeds, shown in Fig. 7. Only radials from the two northernmost HF radars, Hovden and Nyksund, are included due to few drifter positions within the range of the Litløy radar. We find the overall agreement between the two data sources to be good, with a correlation coefficient of 0.72.
As the drifters are rapidly advected northwards after deployment and thus leaving the range of the HF radars, the results discussed in this section are limited to a single data assimilation cycle. The experiments are started on 18 March 2013 at 00:00 UTC. IS4DVAR is run for 24 h, assimilating observations collected during this period and producing modified initial conditions for 18 March 00:00 UTC. Next, a 5-day free simulation is run, initiated from this updated ocean state. In the following discussions, the first day of these simulations is termed "analysis", while the remaining 4-day period is termed "forecast". This is due to the fact that the observations taken within the first 24 h were used during assimilation. It should be noted that the atmospheric forcing used in these experiments is of higher quality than what would have been available in a near real-time setting, which affects the predictability of the ocean forecast.
This procedure has been repeated three times, with different sets of observations to assess the impact the different observation data sets have on the analysis and subsequent forecast. During the first simulation, only CTD hydrography profiles (both temperature and salinity) were assimilated (six profiles just south of the HF coverage area fall within the assimilation window). The observation error standard deviations are set to 0.1 C for temperature and 0.01 for salinity. These are the same values as used by Broquet et al. (2009) and, as argued in their paper, a choice that may enhance the impact of this observation type, which is typically limited in numbers.
In the following, results from this simulation will be denoted "CTD". The second simulation only included HF radar total vectors, and is denoted "HF", while the third simulation assimilated both CTD profiles and HF radar currents, and is denoted "ALL". Additionally, a control simulation, in which no data assimilation was performed, was conducted for the same time period, in the following denoted "CTRL".
During the free simulations, initiated from the analyses produced by IS4DVAR, simulated surface drifters are re- leased in the positions occupied by the real surface drifters from the field campaign. Simulated surface drifters are released every 3 h at the position of the real drifters at that time. The depth of the simulated drifters is set to 65 cm below the sea surface, which corresponds to the average drogue depth of the iSLDMB drifters. In the same manner, numerical floats were released in the control simulation. We validate the results from the analysis period and the forecast period of the simulations separately.
The trajectories of the iSLDMB surface drifters released during the field campaign are compared with their numerical twin from the free simulations following the same methods as used in Röhrs et al. (2012). First, the vector correlation between predicted and observed drifter velocities is evaluated using an analysis similar to the method described in Davis (1985). The trajectories are split into 3 h long segments, and the drift velocities for these segments are then calculated. Defining the vector correlation as the correlation coefficients r for simultaneous pairs of drifter velocities may be calculated. The results are given as a value between 1 and 1, where a value of 1 means perfect correlation in both speed and direction, and a value of 0 means no correlation, while a value of 1 means that the drifters are anti-correlated, i.e. having opposite direction but the same speeds.
The resulting vector correlations between the assimilation simulations and the real drifters are shown in Fig. 8. The quality of the CTRL simulation decreases rapidly and, as the simulation enters into the forecast period, modelled and observed drifter velocities become uncorrelated. The CTD simulation follows the CTRL closely. In fact, the correlation coefficient for the CTD simulation is mostly below the CTRL simulation. As stated above, only six CTD profiles, sampled south of the area where the drifters were deployed, were used. We do not have CTD data co-located with any drifters and, therefore, cannot draw any conclusions regarding the impact of the CTD data. Assimilation of HF radar currents, on the other hand, significantly improves the simulated drifter velocities during both the analysis and the forecast period. Note that adding CTD hydrography to the HF currents further improves the model predictions during the first forecast day, which is in contrast to the detrimental effect of only using CTD observations. A likely explanation is that the addition of hydrographic observations acts as an additional constraint of the vertical density structure, demonstrating the importance of constraining all state variables.
The significance of the improvement has been tested comparing ALL and CTRL using a Wilcoxon rank-sum test. The improvement in current direction is statistically significant, while not so for the speeds. Our interpretation is that the major benefit of assimilating HF radar currents in these simulations is the correction of the current direction, e.g. adjustments in the positions of eddies and the coastal current. This test does not imply that no improvement in speed is obtained (Fig. 8 indicates that the results are better) but that the impact cannot be statistically verified with the limited observational data available.
The impact of data assimilation is also evaluated by comparing observed and predicted trajectories. For this comparison we use the method presented in Liu and Weisberg (2011), in which not only the end points of the observed and modelled trajectories are compared but also the entire history of the drifter trajectories. The normalized cumulative Lagrangian separation is defined as Figure 9. Drifter pathways for two different drifters as observed (green), predicted by CTRL (red) and predicted by ALL (black). The grey tracks shows the pathways of the perturbed initial position floats. The trajectories shown here were released at the start of the analysis.
where d i is the separation distance between observed and modelled trajectory endpoints at time t i after initialization, l oi is the length of the observed trajectory and N is the total number of time steps evaluated. A skill score S is then defined as High skill score means that observed and modelled trajectories agree well throughout the period under evaluation. The skill score is calculated both for the analysis period and for the following 48 h of the forecast period. The results are shown in Table 3. Assimilation do improve predictions of drifter trajectories, although the impact is more limited compared to the drifter velocities. The results show that the skill improves when we consider periods longer than a day. In these cases, data assimilation seems to constrain the ocean circulation in such a way that the predicted trajectories does not stray as far away from the observed paths as they do in a free simulation. As an example of the ability of the model to predict pathways of drifting objects before and after assimilation, two examples of modelled versus observed trajectories are shown in Fig. 9. Note that the drifter in the left panel is released within the area covered by the HF radar antennas, while the drifter in the right panel is released outside of the coverage area.

Comparison with ADCP measurements
We proceed by comparing model results with current measurements from the ADCP rig deployed in our area of interest. A brief comparison of the HF radar total currents with  the ADCP measurements is shown in Fig. 10. Hourly averages of the ADCP observations at 1 m depth have been generated to better match the HF observations, and only observations from periods where observations from both sources are available are included. The ADCP measures greater current speeds than the HF radar, which may be due to the fact that the HF radar currents are spatially averaged as well as temporally. It should also be noted that the dominating current direction is slightly more eastward in the HF measurements and that the current directions measured by the HF radars exhibit greater spread. The results shown in this section are based on 10 sequential data assimilation cycles, where HF radar currents as well as CTD profiles of temperature and salinity are assimilated. Starting from the analysis of the "ALL" experiment described above, the model state at the end of the assimilation window is used as the initial condition for a new assimilation cycle, thus, generating an analysis for the following day. This procedure is then repeated. From each resulting analysis a forecast run is initiated, in the same way as described above. Each forecast thus has a 4-day overlap with the forecast from the previous day.
We compare the results from these forecasts with the . RMSE and bias of transport speed (upper panels) and direction (lower panels) in the water column stretching from 0.5 to 8 m below the surface relative to the ADCP, as a function of forecast day.
To account for discrepancies between the vertical resolution of the model and the observations, we have chosen to focus the comparison on depth-integrated values. To evaluate the depth-integrated flow, we integrate from 0.5 m below the surface and down to 8 m for the upper ADCP, while for the lower ADCP we integrate from 8 m and down to 41 m. Transport magnitude and direction, given by the vertically integrated ADCP currents, are compared with the same quantities calculated from model results. Root mean square error (RMSE) and bias values are calculated as a function of forecast day, and values corresponding to the same forecast day number are binned together From inspection of time series of transport magnitude and direction (not shown), we find that the upper ADCP measures a transport with mostly north-northeasterly heading during the period. The CTRL simulation predicts a northnorthwest-headed current during the period in question. In addition to the discrepancy in direction, the transport magnitude predicted by the CTRL is too weak during the whole simulation in both upper and lower sections. Figures 11 and  12 show mean RMSE and mean bias of the series of simulations along with the means ±1 standard deviation.
In both upper and lower sections it is evident that assimilation causes a reduction in speed RMSE and bias compared with the control run. The change in bias indicates that the predicted current remains more energetic throughout the forecast period. As for direction, the assimilation runs perform better than the control run in the lower section, with reduced RMSE for all forecast days. In the upper section, however, direction predictions seem to be of similar quality in the assimilation runs as in the control run. A possible explanation for this might be that the currents in the near-surface layers are strongly affected by the wind stress. Any changes made to the current in the near-surface layers by the assimilation system may thus have a short life span in the model if the wind stress imposed on the surface yields a different direction. Changes at deeper depths, on the other hand, have a greater chance of being sustained. This advocates for the inclusion of surface forcing in the control vector.

Discussion and concluding remarks
In the presented experiments we have assimilated HF radar and hydrography profiles in both an idealistic and a realistic model. In both cases hydrography profiles were sparse in comparison with the number of HF radar current observations. The realistic experiments were carried out for a specific time period and the results are therefore influenced by local weather conditions and do not necessarily represent the variations in predictability of different flow regimes.
Starting with an idealized test case with synthetic HF total currents and hydrography profiles being assimilated, we find that it is necessary to keep the horizontal error correlation length scales fairly small. The strong currents along the steep topography are associated with sharp fronts and, when large correlation length scales are used, the information provided by the observations may be erroneously spread across these fronts. As the incremental strong constraint 4D-Var is based on repeated iterations of a tangent linear version of the model and its associated adjoint model, the linearization assumption needs to hold throughout the assimilation window. In our case, with strongly nonlinear and energetic currents, this window needs to be quite short. Both the idealized and the realistic experiments give best results when the window is 24 h. We also find that relinearizing the preliminary solution through outer loops is beneficial for the performance of the data assimilation system. These findings are confirmed when similar experiments are performed with a realistic model with real observational data.
Hydrographic observations alone did not yield a significant positive impact on the surface currents, not even in the idealized experiments where more (synthetic) profiles were assimilated, compared to the realistic case. The profiles do, however, have a positive impact on the density field in the model, reducing analysis salinity RMSE from 0.31 in the control run to 0.16 in the experiment where both HF currents and hydrography profiles were assimilated. Also, the temperature RMSE is reduced from 3.0 to 2.2 C. As the ship was moving in the opposite direction of the dominating currents, ship-based observations of temperature and salinity during the forecast period have little value for validation purposes. Comparison with satellite observations of sea surface temperature (not shown here), however, suggests that temperature predictions are improved throughout the forecast when CTD profiles are assimilated.
With assimilation of HF currents, on the other hand, the current forecasts are improved. Gridded fields of surface currents, with temporal resolution of 1 h provide the data assimilation system with improved positioning of eddies and more realistic current speed. In the idealized experiments we even see improvement in the temperature. There are indications that this also holds in the realistic experiments, but independent observations are too sparse to confirm this. Using both HF currents and CTD profiles for assimilation yields the best results. Although CTD profiles did not improve the current forecasts alone, together with the HF radar surface currents they seem to add an additional constraint on the circulation. These results demonstrate the importance of sampling not only the surface but also the subsurface density structure.
We conclude that using HF radar currents in operational data assimilation systems is useful for improving predictions of upper ocean transports, which is highly relevant for oil spill mitigation purposes and search-and-rescue operations. Assimilation of the HF radar radial components directly by utilizing a suitable observation operator holds the potential for further improvement and is the focus of ongoing efforts. Using radials will increase both the number of observations available for assimilation as well as the coverage area. Errors introduced by the algorithms combining radials from two or more antennas will also be eliminated.