Assimilation of HF radar surface currents to optimize forcing in the northwestern Mediterranean Sea

HF radar measurements are used to optimize surface wind forcing and baroclinic open boundary condition forcing in order to constrain model coastal surface currents. This method is applied to a northwestern Mediterranean (NWM) regional primitive equation model configuration. A new radar data set, provided by two radars deployed in the Toulon area (France), is used. To our knowledge, this is the first time that radar measurements of the NWM Sea are assimilated into a circulation model. Special attention has been paid to the improvement of the model coastal current in terms of speed and position. The data assimilation method uses an ensemble Kalman smoother to optimize forcing in order to improve the model trajectory. Twin experiments are initially performed to evaluate the method skills. Real measurements are then fed into the circulation model and significant improvements to the modeled surface currents, when compared to observations, are obtained.


Introduction
The ocean coastal dynamics of the northwestern Mediterranean (NWM) Sea has been studied over the last decades using observations (e.g., Millot, 1999;Petrenko, 2003;Sammari et al., 1995;Allou et al., 2010), modeling (Echevin et al., 2003;Casella et al., 2011;Ourmières et al., 2011;Pairaud et al., 2011), or a combination of both (Schaeffer et al., 2011b;Guihou et al., 2013). The NWM is a micro-tidal sea (e.g., Alberola et al., 1995) and its dynamic is mainly governed by wind forcing and by the Northern Current (NC), which is the northern branch of the general circulation of the NWM Sea (Millot, 1999). The NC is a quasi-permanent density current that originates in the eastern part of the NWM and forms a large-scale loop. Bathymetric constraints operate on the northern, western and southwestern sides of its path. This geostrophic current can experience a strong variability in this particular region, due to the steep topography and the irregular coastline, as well as the atmospheric forcing. The NC varies in strength and width, its meandering can generate eddies, and its separation in the vicinity of Toulon ( Fig. 1) controls water fluxes at the entrance of the Gulf of Lion (GoL). Several of the above studies were carried out in the GoL or off Nice. However, the dynamics of the NC in the intermediate area, typically the region off Toulon, is poorly documented. Since May 2010, an HF radar system has been deployed in this region in an attempt to fill this void, allowing the gathering of continuous and synoptic observations of surface currents up to 100 km off the coastline with typical temporal and spatial resolutions of 1 h and 3 km, respectively.
Scientific interest in the assimilation of HF radar measurements in hydrodynamical models has been growing over the last few years. Various assimilation methods have been developed based on nudging (Lewis et al., 1998), Kalman filtering using parametrized (Paduan and Shulman, 2004;Shulman and Paduan, 2009) or ensemble (Breivik and Saetra, 2001;Oke et al., 2002b;Barth et al., 2008a) error covariances, and the application of adjoint-based methods (Kurapov et al., 2003;Hoteit et al., 2009).
The method applied herein consists in optimizing the wind-forcing and lateral boundary-forcing components of a . Red points mark radar stations PEY and BEN. The black contour lines show areas where HF radar data are available more than 20 % over the assimilation period for PEY (full) and BEN (dashed) radars and is called the observation domain (OD). The black box sketches the pseudo-observation domain (POD). The meridional magenta line crossing the coastal current will be used for Hovmöller diagrams. Surface currents of free simulation averaged over the assimilation period are shown. Color bars range from 0 (white) to 0.5 (red) m s −1 . The 100, 1000 and 2000 m isobaths are drawn (black thin lines).
numerical regional primitive equation model in order to minimize differences between the model and the radar-derived surface currents.
Wind optimization is an important issue as, in coastal zones, surface wind errors can explain a large part of the model error regarding surface currents (He et al., 2004;Barth et al., 2008b). Such optimization can be achieved by correcting the surface wind field (Barth et al., 2011) or the parameters governing the wind-sea momentum transfer, e.g., the drag coefficient (Ten Brummelhuis et al., 1993;Peng et al., 2012).
Similarly, open boundary condition (OBC) optimization can be achieved by correcting the OBC values themselves or by correcting the OBC's formulation and discretization. This last approach has been addressed several times (e.g., Orlanski, 1976;Flather, 1976;Marchesiello et al., 2001;Blayo and Debreu, 2005). It consists in a numerical treatment of the oceanic boundaries allowing the removal of the perturbations generated inside the model domain without deteriorating the inner model solution (Røed and Cooper, 1986). Research in this field is still ongoing. Instead of pursuing this approach further, we shall address the first approach, i.e., the optimization of time-dependent OBC values. This has been achieved previously in barotropic models (e.g., Seiler, 1993;Ten Brummelhuis et al., 1993;Gunson and Malanotte-Rizzoli, 1996;Bogden et al., 1996;Taillandier et al., 2004) and more rarely in baroclinic models (e.g., Shulman et al., 1999;Taillandier et al., 2009).
Variational methods (e.g., Le Dimet and Talagrand, 1986) are often preferred to tackle forcing optimization problems. The adjoint model technique is the one most often used, but depending on the model it can sometimes be difficult to formulate (Shulman et al., 1999). Another way of solving optimization problems is by the application of a Kalman filter (e.g., Gelb, 1974). However, in its original form, this technique has a prohibitively high computational burden to be used in geophysical applications. The development of reduced rank methods such as ensemble methods (Evensen, 1994) allows for the reduction of this computational cost because only a limited number of model states (O(100)) needs to be stored in order to obtain an accurate estimate of the error covariance matrices. Moreover, ensemble methods are efficient for dealing with error propagation in non-linear models.
Recently, the Ensemble Perturbation Smoother (EnPS) has been developed to optimize stationary tidal boundary conditions of a hydrodynamical model of the German Bight (Barth et al., 2010), a very shallow area (maximum depth of about 50 m) where winds and tides are the main sources of variability. The EnPS has also been used successfully to optimize surface wind estimates (Barth et al., 2011) in the same area, but the resulting model surface current was not directly evaluated. This beckons the following questions: firstly, whether optimizing wind can improve model surface currents, and secondly, whether it is possible to optimize non-stationary baroclinic OBC by applying the EnPS with the aim of improving model surface currents in a highly dynamic area. These issues are addressed by using the EnPS in a regional baroclinic model of the NWM Sea and by using the HF radarderived surface currents measured off Toulon.
Section 2 describes the hydrodynamical model configuration and the experimental data set. The data assimilation method is explained in Sect. 3. Results of wind optimization are shown in Sects. 4 and 5, showing the results of OBC optimization. Finally, results are discussed in the conclusion.

Model
The simulations described below were obtained using the Regional Ocean Modeling System (ROMS, McWilliams, 2003, 2005), which solves the ocean primitive equations, featuring a free surface formulation and hydrostatic approximation. It uses terrain-following coordinates in the vertical dimension and curvilinear coordinates in the horizontal dimension. ROMS incorporates advanced features and high-order numerics, allowing efficient and robust resolution of mesoscale dynamics in the oceanic and coastal domains. ROMS has been widely applied to different oceanic regions for a variety of ocean dynamics studies (Penven et al., 2001;Capet et al., 2008;Casella et al., 2011). The model domain (MD) as depicted in Fig. 1  Two active open boundaries connecting the MD with the NWM basin are incorporated. The OBC treatment is a well-known ill-posed problem (e.g., Chapman, 1985;Marchesiello et al., 2001;Blayo and Debreu, 2005), and the goal here is to implement the optimal conditions in order to allow inflow from the open sea, to conserve mass and energy and to reduce any wave reflection phenomena. A radiation-relaxation condition for 3-D momentum and tracers is combined with a Flather condition for the 2-D momentum (Marchesiello et al., 2001), and a sponge layer is added to avoid spurious reflections. The typical relaxation coefficients have been adjusted to 1 and 10 days (for tracers and baroclinic velocities, respectively) for inflow and to 360 days for outflow.
The data for the relaxation boundary conditions were derived from the PSY2V4R2 operational configuration developed by MERCATOR-OCEAN (http://www.mercator-ocean.fr), based on the NEMO (Nucleus for European Modelling of the Ocean) model. The PSY2V4R2 configuration features a horizontal resolution of 1/12 • with 50 vertical z levels, and includes a sequential data assimilation system (SAM2V1) based on the Kalman filter. Satellite sea surface temperature (SST) and sea level anomaly (SLA) are assimilated, in addition to available in situ temperature and salinity profiles. Daily averages of the velocity, mass fields and surface height at the MD open boundaries are used. The forcing is made off-line.
At the air-sea interface, the oceanic model is forced by the ALADIN operational regional model from Météo-France. ALADIN runs at a resolution of 1/10 • and features data assimilation and state-of-the-art atmospheric physics (Fischer et al., 2005). The importance of such high spatial and temporal resolution in terms of oceanic response has been widely studied and demonstrated over the last few years, particularly to achieve a more realistic representation of the oceanic mesoscale circulation (Langlais et al., 2009;Schaeffer et al., 2011a).
The initial condition of the simulation is based on the daily averaged baroclinic and barotropic velocities, temperature, salinity and surface elevation derived from PSY2V4R2 for 1 July 2011. A 4 month spin-up period is considered to be adequate for reaching a state of equilibrium. The model output corresponds to 3-hour averaged fields.
The average surface current for the period 15 November-15 December 2012 is shown in Fig. 1. The model simulations give a maximum mean current of 0.5 m s −1 , in accordance with the observations.

Observations
Since May 2010, HF radars have been monitoring surface currents off Toulon (Fig. 1). The radars employed are WERAs (Wellen Radar, Gurgel et al., 1999) and operate at 16.15 MHz with 50 kHz bandwidth. A radar system is installed at Fort Peyras in Cape Sicié (hereafter PEY) and operates in a monostatic configuration, i.e., transmitter and receiver are very close together (∼ 100 m). The second system operates in a bistatic configuration. The transmitting antennas are those of the PEY radar, whereas the receiving antennas consist of a linear array of eight antennas installed near Cape Bénat (hereafter BEN), 40 km to the east of PEY. Integration time, i.e., the time over which the temporal signals are processed, is one hour. Velocity accuracy is 4 cm s −1 , and the spatial sampling is 3 km in range. An azimuthal discretization of 2 • was achieved through direction finding using the MUSIC processing technique (Schmidt, 1986). This technique is now in common use in the HF radar community (Lipa et al., 2006;Sentchev et al., 2013). The radial component of the current is measured along the look direction of the radar. With the given radar frequency, the radial velocity is theoretically measured at a depth of 74 cm (Stewart and Joy, 1974).
The MUSIC technique involves spatial holes, temporal gaps and false alarms (outliers) on radial velocity maps. For the present study, a series over a 30 day time period was used at every radar cell in order to remove outliers. Outlier removal was performed using the histogram of the temporal gradients of the current. Spatial holes and temporal gaps were not filled, in order to avoid the addition of spu- rious and correlated information in radial velocities, which can result from interpolations. Only the radial component maps were taken into consideration, and not the current vectors that can be deduced by the geometric combination of the two sets of radar measurements. The first reason for this choice is to prevent the introduction of the geometric dilution of precision, which is inherent in the vector mapping process (Chapman et al., 1997). The second reason is to dispose of the most extended observation area and to avoid wasting valid observations. HF radar measurements from PEY and BEN (y pey , y ben ) gathered during the period 15 November 2011 00:00 UTC to 15 December 2011 23:00 UTC were taken into consideration for this study. A spatial mask rejecting the input from radar cells having a temporal coverage below 20 % over the period was applied to the data. The 20 % temporal coverage level is depicted in Fig. 1 and is hereafter referred to as the observed domain (OD). When measurements from both stations are available, variability in the total coverage is mainly due to the degradation of PEY performance owing to high wind speed.
Twenty-eight temperature and salinity profiles from conductivity-temperature-depth (CTD) measurements were also taken into consideration for the analysis of assimilation results. CTD measurements were carried out in the area within the geographical grid defined by coordinates 5.8-7.35 • E and 42.6-43.3 • N, over the period 11-14 December 2011, during a cruise undertaken in the framework of the European TOSCA (Tracking Oil Spill and Coastal Awareness network, http://www.tosca.med.eu) program. Note that this data set is independent data that can be used for validation of results.

Model validation
The Hovmöller diagram (time versus latitude) of both modelderived and observed radial velocities is depicted for the PEY station (Fig. 2). The radial velocity is the component of the total surface velocity of the current U in the radial direction. By convention its sign is positive (negative) for velocity going toward (away from) the station. High posi-tive values in the PEY radial component are a distinct signature of the coastal current (Marmain et al., 2011). The model clearly underestimates the intensity of the current radial components, and the position of the maximum seems slightly shifted off-shore. However, the temporal variability is well simulated: the change in sign at the beginning of December, as well as the high-frequency modulations (near-inertial oscillations), and the maximum intensity observed from PEY around 23 November are reproduced well by the model.
The validity of the vertical structure of the model was evaluated using CTD measurements. Model temperature and salinity exhibit low departures from observations. According to the definition of the root mean square error (RMSE), Er and the bias B given in Sect. 3, Er = 0.780 • C and B = 0.002 • C for temperature and Er = 0.084 and B = −0.054 for salinity. Significant correlations appear (0.79 and 0.86 respectively). The main errors come from the inaccurate representation of the thermocline by the model. However, densities are in good agreement.

Data assimilation methodology
HF radar-derived velocities have already been used to correct M2 tidal boundary conditions (Barth et al., 2010) and surface wind fields (Barth et al., 2011) using the so-called EnPS, a fixed-interval sequential smoother (Cosme et al., 2012). In this study, we investigate, firstly, whether optimizing wind can improve model surface currents, and secondly, whether it is possible to optimize non-stationary baroclinic OBC by applying the EnPS with the aim of improving model surface currents in a highly dynamic area. Uncertainties in these atmospheric and oceanic fields are represented by an ensemble of wind forcing or an ensemble of OBC forcing based on perturbations of the original (or background) fields (x b ).

Ensemble setup
The perturbation method is detailed in Barth et al. (2011). Briefly, the spatial structure of the forcing perturbations is obtained by using a Fourier decomposition of the forcing field. The Fourier coefficients correspond to the spatial field. For each temporal frequency subjected to perturbation, a complex random time series is constructed with a temporal correlation scale according to the perturbed frequency, with nil mean and unit variance (Evensen, 1994). Such a series is then multiplied by the corresponding Fourier coefficient and the resulting products are summed. The real part of this quantity corresponds to the forcing perturbation x l p , with l as the member index. x l p is added to x b , resulting in With this method, forcing errors are assumed to have similar spatial and temporal scales to the original, non-perturbed forcing. By using this approach, it is assumed that the uncertainty in the forcing field is mainly caused by phase and amplitude errors of the underlying processes rather than errors in the variability. Perturbations are estimated using the time variability of the forcing field. x l p is weighted by the α coefficient in order to achieve realistic perturbations. This approach ensures that a spatial structure is multiplied by a compatible temporal scale. The forcing ensemble must verify x b = x l , with . as the ensemble average operator. An ensemble of 100 perturbations for zonal and meridional components of the wind vector W was generated over the MD and over a period extending from 1 November to 31 December 2011, taking into consideration temporal scales ranging from 3 h to 30 days. The minimum period is the temporal resolution of the original forcing and the maximum is set to 30 days, assuming that higher periods do not contain uncertainties.
An ensemble of 100 OBC perturbations was generated, taking into consideration temporal scales ranging from 1 to 30 days, over the same time period, for temperature, salinity, surface elevation, zonal velocities and meridional velocities (respectively, T obc , S obc , η obc , U obc , V obc ) at the eastern boundary of the MD. Only this boundary was taken into consideration because it is located upstream of the observation area along the NC flow. Perturbations were introduced in the northern part of the OBC (42.5-44 • N) and over the sub-surface ocean layer (−200-0 m) as the prognostic variables at higher depths do not exhibit any significant correlation with the surface current in the OD. In fact, such a restriction on the perturbation domain can be seen as a spatial localization (Hamill et al., 2001;Houtekamer and Mitchell, 2001). Therefore, we consider that errors in the model solution originate only from the perturbed area.
For each forcing perturbation, the model was run using the same initial condition as the free simulation. A 15 day ensemble spin-up was selected to attain a satisfactory variance between surface current members. The ensemble simulations for the time period 15 November-15 December 2011 were taken into consideration.

A general description of the EnPS (Ensemble Perturbation
Smoother) is given in Appendix A. For more information, the reader is referred to Barth et al. (2010Barth et al. ( , 2011. Since the EnPS is a non-sequential method, all available observations over the period of analysis are used to compute the optimal state. The time dimension is consequently embedded in the observation and optimization vectors. The model trajectory is not estimated directly. Rather, the optimal forcing trajectory is derived and then run in the model. Furthermore, the advantage of estimating the forcing first is that the final model trajectory satisfies (per construction) the model equations (and thus all non-linear balances between model variables), which is not the case if the model trajectory is analyzed directly, and correcting the forcing instead of the model state also avoids unrealistic transients (Malanotte-Rizzoli et al., 1989).
The optimized forcing x a is computed using the Kalman filter analysis step with a non-linear observation operator: where K is the Kalman gain that varies according to ensemble and observation error covariance matrices (Appendix A). The non-linear observation operator h(.) is applied to the forcing x and h(x) represents the surface currents within the integration period. The analyzed model trajectory is obtained in a final step, by performing a new run of the model using the optimized forcing field.
The analysis x a is unbiased if both background estimate x b and observations y o are unbiased. If the observation operator is non-linear, it can be verified that the analysis is unbiased using Eq. (3) (Barth et al., 2011): The degree of linearity of the observation operator was verified using Eq. (3). The overall root mean square difference between both sides of Eq. (3) was computed for wind and OBC ensembles. A value of 0.03 m s −1 was found, which is small in comparison with other errors and bias. However, the standard analysis methods provide often useful results even if the background is biased. The observation vector y o is associated with an observation error covariance matrix R. For the sake of simplicity, R is assumed diagonal and it is also assumed that R is the sum of the measurement error covariance R obs I and the representativity error covariance S 2 obs I where I is the identity matrix. R obs represents the measurement errors that arise from noise effects imparted by data processing as well as by the radar itself. S obs is introduced to characterize processes and scales present in the observation but not resolved by the model. The representativity error could be extended to include processes that cannot be corrected using these forcings only (Barth et al., 2011). To increase the effective dimension of the error subspace, to discard distant observations involving unrealistic correlations and to improve computational capabilities, the Kalman filter analysis step can be performed in reduced spatial domains, zone by zone, which is known as spatial localization (Hamill et al., 2001;Houtekamer and Mitchell, 2001). A similar approach can be taken in the temporal dimension. Typically, the model integration period is divided into a series of shorter, non-overlapped, intervals called optimization windows with a duration of win ana . Considering an observation window of duration win obs , such as win ana = win obs , forcing fields within the time range [t; +win ana ] are updated using all available observations in such observation windows, with t as the beginning of the optimization window (Barth et al., 2010(Barth et al., , 2011. Here, the definition of the observation window was modified to take into account backward (win − obs ) or forward (win + obs ) observations within the time range t − win − obs ; t + win + obs (Fig. 3), considering the delay of the oceanic response to the lateral and surface forcing that may result in non-instantaneous correlations between forcing and surface currents (e.g., Seiler, 1993). This way, the standard EnPS is retrieved with win + obs = win ana and win − obs = 0.

Experimental setup
The assimilation methodology was first applied to a set of radar pseudo-observations in identical twin experiments (TE) to assess the method, and then applied to experimental radar observations. In TE, forcing and model control states (x true , h(x true )) are extracted from the ensemble members. These control states are used to generate a set of pseudo-radial velocity fields by adding a Gaussian noise. The free run corresponds to running the model using the original forcing. The assimilation of these pseudo-observations results in an analyzed forcing that can be used in the model to simulate the resulting analyzed model surface current. The true state is known and the skill of the assimilation scheme can be quantified.
Pseudo-observation values are given for a geographical grid defined by coordinates 5-7 • E and 42.2-43.2 • N, the pseudo observation domain (POD), and are defined as the projection of model surface velocities in the radial directions at each model node. The measurement noise is assumed to be nil.
The experimental observations consist in the radial velocity measurements derived from the two PEY and BEN radars (Sect. 2.2) with R obs is 0.04 m s −1 , corresponding to the velocity accuracy.

Diagnostics
The skill of the assimilation method is evaluated by comparing the control state (control) to the free (b) or analyzed (a) states in terms of RMSE (e.g., Oke et al., 2002a). In the case of the twin experiments, the control state is the true one, true. In the real case experiments, a cross-validationlike technique is used (e.g., Brankart and Brasseur, 1996): the control state is a random set of observations, o, discarded at each time step from the initial set of observations.
Let n be the spatial index of dimension N, and k the temporal index of dimension K. The RMSE of a variable is defined by When the summation is applied to both spatial and temporal dimensions of , Er corresponds to an overall scalar measurement of the assimilation performance. If the sum is computed in the spatial dimension only, the error is evaluated during the experiment time evolution (Er N (k)), and when the average is calculated in the temporal dimension, the performance is evaluated in terms of spatial structure (Er K (n)).
The bias B can also be computed to determine the mean difference between model solution and control: The skill score "Sc" is an additional metric that can measure the performance of the data assimilation (e.g., Oke et al., 2002a), and is defined by A positive score "Sc" indicates that the assimilation is efficient, while a negative score indicates a degradation. In particular, Sc = 1 is obtained for a perfect analyzed field, and Sc = 0 reflects a nil impact on average. Eventually, = Er b K (n) − Er a K (n) will also be computed.

Results of a twin experiment
The first sensitivity study concerns the wind forcing. Wind perturbations are introduced in such a way that the resulting wind ensemble represents the wind uncertainties. Insufficient wind perturbations may induce only a small variance of the model surface current ensemble and the model control space will be misrepresented by this ensemble of model states. On the contrary, excessive wind perturbations can generate unrealistic patterns in the wind ensemble or in the corresponding model ensemble. We found that, for the present model setting, an error in wind forcing corresponding to 50 % of its temporal variability (α = 0.5) is the most efficient in terms of the overall RMSE between control and analyzed states for wind and model solutions during the assimilation period.
The representativity error, S obs , is also investigated. We look for the value of S obs that provides the most accurate optimized wind and model solutions for a given optimization and observation window. A high value of S obs can give insufficient weighting to the observations relative to the background state, while a low value of S obs can produce unrealistic optimized forcings. Therefore, the analysis of the forcing and the model solution must be investigated concurrently to find the optimal compromise. For the present experiment, S obs was set at 1 m s −1 .
The assimilation process also requires the selection of the duration of the temporal optimization window, win ana , used to localize the analysis. This duration should minimize the discontinuities between two consecutive windows. Typically, the amplitude of these discontinuities must be on the same order as the temporal variability of the model, and again a compromise must be found to reduce the error of optimized forcing and the resulting model solution. 1 day was found to be adequate for win ana .
Traditionally, observations are considered within the optimization window. However, depending on the application, the covariances between observations and ensemble simulations at a given time may be low due to a weak causality between forcing and model solution at the observation points. At such a stage we suggest taking into consideration a temporal shift between the observation and optimization windows. Given a time reference that corresponds to the beginning of win ana , the observations before (after) this time are included in win − obs (win + obs ). The values win − obs = 0 day and win + obs = 3 days provide the best results in terms of the overall Er for wind and surface currents. The experiment performed with this data set is hereafter referred to as WindTE. Figure 4a shows the improvement in wind field estimation using idealized radar measurements in terms of the time evolution of Er N (k). The assimilation is effective over the total period except at some dates (e.g., 30 November-2 December 2011), representing only 8.6 % of the overall duration. These dates mainly correspond to weak covariances between surface currents on the POD and wind over the MD, corresponding to the presence of mesoscale circulation patterns that are not wind induced.
The resulting velocity field in the radial component shows an improvement during the assimilation period after a 3-day adjustment period (Fig. 4b). The overall Sc of radial velocities is 0.54. Oscillations can be observed to have a period of 17 h, which is close to the inertial period at the considered latitudes, but cannot be directly correlated with the wind solutions (Fig. 4a). The presence of such oscillations indicates that free and control surface current near-inertial oscillations are not in phase and that the wind optimization has only a limited impact on the near-inertial oscillations.  The spatial patterns of Er K (n) for the forcing wind are depicted in Fig. 5, superimposed on the negative values of to highlight the regions where the optimized wind is degraded compared to the original wind. Over most of the MD, assimilation of surface currents improves the wind estimation. The area in the eastern part of the MD, where wind estimation is degraded, corresponds to a region where original and control winds are in good agreement before assimilation.
The model response to the analyzed wind forcing in terms of Er of surface velocity field is shown in Fig. 6. The surface current vector field resulting from analysis is generally more accurate than the free run. Maximum improvement is observed at the entrance of the GoL and along the shelf break.
Along the eastern open boundary the assimilated solution is less accurate, corresponding to the degradation of the analyzed wind forcing (Fig. 5).
Hovmöller diagrams of radial velocities along a southnorth transect at 6.2 • N across the POD are shown in Fig. 7 for the PEY radar station. The efficiency of the optimized wind forcing on the radial velocity model response is visible in the high intensity patterns. For the positive radial velocities (red), the assimilation reduces the overestimation of the free run in the period 20-24 November. Later in December, the free run clearly underestimates the radial velocities, and again the optimized wind intensifies these features, for outward components (4-10 December) and inward components (10-16 December). Moreover, the use of the optimized wind preserves the near-inertial oscillations. Results for the BEN radar station (not shown) are similar but weaker.

Results of a real case experiment
The WindReal experiment is performed to correct the same free model run with real radar observations, from the two PEY and BEN radar stations (Fig. 1). Observations are averaged over consecutive 3 h windows to improve computational capabilities.
The best compromise in terms of lowest overall Er and Sc is obtained for α = 0.3, S obs = 1 m s −1 , win ana = 1 day, win − obs = 0 day and win + obs = 3 days, according to the findings from sensitivity tests performed on these parameters.
For this experiment, contrary to the TE where the true state is known in terms of wind and model solution, no ground truth data are available to validate the optimized wind forcing. Meteorological stations in the region did supply wind observations effectively, but all were assimilated into the AL-ADIN meteorological forcing, and no additional independent wind observations were available. To ensure realistically optimized winds, even if this study does not aim at improving wind fields, the original and optimized winds have been compared to in situ winds from four land stations located near Toulon and from two moorings located in the GoL and off  The time evolution of Er N (k) computed from the radial surface current components measured by the PEY and BEN radars is plotted in Fig. 8.
Current estimation by using wind optimization is better for the PEY station than for the BEN station, with an Sc of 0.36 and 0.13, respectively (Fig. 8b, c, and Table 1). The two radar stations have different viewing geometries, in particular regarding the NC vein, which is the most dynamically active region where errors are the highest (Fig. 9a, d). While the PEY radar station is well placed to identify the westward coastal current velocity corresponding to the largest component of the NC, the BEN station mainly measures its mean-dering and oscillations, as the NC flows perpendicularly to the look direction of the BEN radar. Table 1 shows radial velocities Er, B and Sc for PEY, BEN and both radar stations. When both stations are considered together, Er reaches values of up to 0.15 m s −1 (Table 1), and the overall bias is halved.
The spatial patterns of Er K (n) are limited to the radar regions ( Fig. 9) to be compared with the observations. Free radial velocities from the PEY station show higher discrepancies with the observations (Fig. 9a) than free radial velocities from the BEN station (Fig. 9d). The assimilation performs well in most of the domain, especially north of 42.7 • N. The best improvement (red areas of the figure) corresponds to the NC vein, the most energetic feature of the region where errors are the highest.
Hovmöller diagrams from free (Fig. 2a) and analyzed (Fig. 11a) radial velocities from the PEY station show that optimized wind forcing acts on the NC vein intensity. During the period 15-26 November 2011, the intensity of the NC vein is strengthened, better approximating the observations (Fig. 2b), and is weakened during the period 4-9 December 2011. A southward shift of the NC vein was also found during the period 10-16 November 2011.
Data assimilation performance reduces when moving southward and degradation regions may appear. This may be due to the scarcity of observations when moving away from the radar stations.
The norm of vector speed difference U b −U a between free currents U b and analyzed currents U a averaged over the assimilation period and at the surface (Fig. 10a) and integrated over the first 200 m below the surface (Fig. 10c) shows that the corrections are primarily located on the surface in the NC vein all along the French coast and over the GoL shelf break and can occur over the whole MD. However, they vanish with depth.
Once we have verified that the wind optimization obtained through the assimilation of surface velocities does  not degrade the model solution in the observation area, the question of its impact off the OD, and in the whole water column, in terms of hydrological variables, arises. Temperature and salinity from CTD measurements and free/analyzed simulations are compared. The wind optimization method employed only has a weak effect on the temperature (B = 0.017 • , Er = 0.818 • ) and salinity (B = −0.049, Er = 0.081) profiles. No significant changes in correlation be-tween observed and free or WindReal analyzed temperature and salinity are noted. This result is not surprising, as the primary wind effect on the ocean dynamics is surface current response, and only for longer timescales should the turbulent processes transfer the surface instant impact to deeper layers. However, assimilation of radial velocities may not enable the corresponding wind forcing to impact the whole water column or any hydrological variables.

Results of a twin experiment
An attempt to correct the OBC values using the same method as for the wind correction is evaluated first on the TE. These experiments are performed with 1 day-smoothed atmospheric forcing to weaken the influence of the highfrequency wind variability on surface currents and focus on the OBC influence only.
The method requires the tuning of the same parameters as for wind optimization: α, S obs , win ana , win + obs , win − obs . The reference control state is a member of the ensemble with α = 1. After several numerical tests, the best compromise in terms of the overall Er and Sc was found with α = 1.5, S obs = 2 m s −1 , win ana = 1 day, win − obs = 0 day and win + obs = 15 days. The experiment performed with this data set is hereafter referred to as ObcTE.
It should be noted that the choice win +,− obs is crucial. As the coastal circulation flows westward, win + obs allows the temporal localization to take into account the delay necessary for propagating OBC information up to the observation area. This delay is about 15 days for an advection propagation of 25 cm s −1 (e.g., Fig. 1). The nil value obtained for win − obs denotes that no valuable information propagates from the observation area to the eastern OBC.
The assimilation method leads to a significant improvement in optimized OBC variables U obc and η obc . The method   has only a limited effect on T obc and S obc , while V obc is degraded (Table 2). This suggests that zonal velocity and surface elevation at the eastern boundary are the key parameters controlling model surface currents in the POD. The first parameter imposes the fluxes at the entrance of the MD, whereas the second parameter governs the geostrophic component of the NC, which is the main large-scale feature of the region (e.g., Birol et al., 2010). The degradation of the meridional velocity reflects the weak correlation that exists between surface currents in the POD and V obc . The spatial pattern of Er a,b K (n) at the eastern boundary of the model shows that the maximum improvement for U obc and η obc occurs north of 43.5 • N (not shown).
In the POD, the optimized OBC values improve the model radial velocity (Fig. 12a). From 15 November to 18 November 2011, no effect on the surface currents is apparent. This corresponds to the time necessary for the OBC information to reach the POD. From 25 November 2011, Er a N (k) is smaller than Er b N (k), with differences of up to 1 cm s −1 , except at the end of the time period. The lower efficiency of assimilation on surface currents from 14 December 2011 is a direct consequence of the degradation of U obc . The overall Sc of radial velocities is 0.22. Similar comments can be made for the vectorial surface currents (Fig. 12b), except that no degradation occurs at the end of the time period.
Despite several regions of degradation, model radial surface currents are globally improved, as well as surface temperature and salinity (not shown).

Results of a real case experiment
The same free run and the same real radar observations as in WindReal are used to correct OBC values. In order to set the optimal value for α, the nearest observations were used to get an estimate of the forcing error relative to its temporal variability at the eastern boundary. α ranged from 0.5 to 1. In this study, the lowest overall Er were obtained for α = 1. The value S obs = 2 m s −1 was the best compromise for the representativity error. The other parameters are identical to those of the ObcTE experiment. The experiment performed with this data set is hereafter referred to as ObcReal. Figure 13 shows the temporal evolution of Er N (k) for radial surface currents, for each individual radar site, and for both. When both stations are considered, the analyzed currents better approximate observations, with differences between Er a N (k) and Er b N (k) of up to 7 cm s −1 . The adjustment delay at the beginning of the period is needed for the optimized information to be propagated from the eastern boundary to the OD. Statistical analysis reveals that the absolute value of the overall B is reduced by 2.3 cm s −1 for the analyzed surface currents compared to the free surface currents ( Table 1). As in WindReal, greater improvements are obtained for the PEY station (Fig. 13b) and the impact of data assimilation is small on BEN radial velocities (Fig. 13c).
Er K (n) maps for PEY and BEN are shown in Fig. 9c, f. Modeled radar currents are mainly improved north of 42.8 • N in the NC vein. For PEY, no significant degradation is observed. On the contrary, for BEN, areas of degradation up to 1.6 cm s −1 are observed around 42.6 • N, which corresponds to the southern edge of the NC. The analysis of vectorial surface currents (not shown here) shows that, using the optimized OBC, the NC vein is more marked due to higher surface velocities along the whole NC pathway.
Hovmöller diagrams of free (Fig. 2a) and analyzed (Fig. 11b) radial velocities from the PEY station show that the use of optimized OBC values has a similar effect on radial velocities to the use of optimized wind fields but to a lesser extent.
The structure of corrections is investigated by computing the norm of U b − U a averaged over the assimilation period, at the surface (Fig. 10b) and integrated over the first 200 m below the surface (Fig. 10d). It appears that corrections are primarily located in the NC vein all along the French coast and over the GoL shelf break and affect a large part of the water column. This result proves clearly that the NC strongly depends on its upstream dynamic.
The effect of OBC optimization on the water column is also studied using CTD measurements. Correlations between observed and free or OBCreal analyzed temperature and salinity remain unchanged. A bias difference between ObcReal analyzed and free runs of 0.12 • C is observed in temperature values, which indicates that the effect of assimilation on the water column is more important than in Win-dReal. A slight degradation is obtained for salinity values, with a bias difference of −0.006.

Conclusions
For the first time, HF radar measurements in the NWM Sea are assimilated into a circulation model. The region is characterized by a steep topography and complex coastline, and its ocean dynamics includes a geostrophic coastal current associated with strong mesoscale activity. The EnPS is used to correct two different sources of error acting on the model surface current: high-resolution wind forcing and baroclinic open boundary conditions.
The performance of the assimilation method can be assessed by evaluating either the optimized forcing or the resulting model response. Since no independent wind measurements were available, the assessment of the wind optimization could be done only on the TEs, which clearly demonstrate the strong non-local impact of this method. A limited area of HF radar observations can be effective in large-scale wind correction, and it can reasonably be inferred that the resulting model response will be improved over the whole area. This method also improves surface current simulations through OBC optimization. This result is non-trivial because the OBC optimized area is located far from the observation area (150 km). The perturbation of OBC is done in the NC vein area and takes into account uncertainties regarding width, depth and speed of the NC vein since all variables are perturbed.
The results show that the correction of the wind forcing has a stronger impact on the resulting simulated surface current than the OBC forcing correction. WindReal Er is 1 cm s −1 less than the ObcReal one, and in both cases, bias is reduced. Wind optimization effects are felt in surface currents over short time intervals contrarily to OBC optimization, whose improved values are propagated by model equations over time and space. In terms of temperature and salinity along the water column, mixed results are achieved. No significant improvements are obtained and areas of degradations can occur in temperature values for WindReal or salinity values for ObcReal.
The assimilation of HF radar-derived radial velocities acts on model surface currents over the whole of the MD. The main corrections are depicted by red values in Fig. 10, which represent the norm of U b − U a averaged over the assimilation period. These corrections are primarily located in the NC vein all along the French coast and over the GoL shelf break. WindReal corrections are less significant off Nice but occur over the whole of the MD. When the same vector speed difference is applied to the 200 m-depth averaged currents, the correction introduced by ObcReal can affect model currents of the deeper water masses, whereas the correction introduced by WindReal mainly concerns model surface currents.
These results inspire confidence in the implementation of the method on OBC forcing optimization. When dealing with a domain with large open boundaries that strongly influences internal ocean dynamics, a small inconsistency from the outside will certainly trigger unrealistic features in the domain of interest. The availability of local monitoring systems in the area of interest could be used with this method, to correct the OBC input information, in terms of intensity or positioning of incoming/outcoming currents. In addition, the monitoring of the region off Toulon can be of value for the comprehension of the NWM sea circulation, in particular when considering the NC variability, which is still poorly documented. The recent development of HF radars on the Mediterranean coast through national and international projects should be extended to enhance the coastal monitoring and the assimilation of HF radar data.
The methodology was successfully applied and differences between observed and modeled surface currents were significantly reduced. However, further improvements can be envisaged by taking more aspects into consideration. Wind and OBC optimization have been tested independently and the impact on the model by simultaneous forcing optimizations should be considered in future work, which would require the use of a more realistic eddy-resolving model. Only the eastern boundary was addressed here, and the complementary treatment of the southern boundary may also improve the overall result. Furthermore, the regional model used in this study is quite coarse and may contain errors not directly linked to the forcing: for instance, the bathymetry and the coastline contour may not be well resolved, or the parametrization of the OBC treatment or the turbulent schemes may not be accurate enough when considering the assimilation of high-resolution surface currents derived from HF radar. The use of a high-resolution regional model configuration with an advanced physical parametrization (Ourmières et al., 2011;Guihou et al., 2013) could also lead to significant improvements.

Appendix A: Ensemble Perturbation Smoother
EnPS was proposed by Barth et al. (2010) to optimize tidal boundary conditions and was successfully applied by Barth et al. (2011) to optimize surface wind forcing. Contrary to the classical sequential method, which intends to derive the best estimation of the model state, this method aims at optimizing forcing values. These forcing fields are then run once again in the model to get the final model solution. Thus, the model trajectory is corrected instead of the model state as in the ensemble smoother (Van Leeuwen, 2001), the 4D-LETKF (Hunt et al., 2004(Hunt et al., , 2007 or the AEnKF (Sakov et al., 2010). This method also ensures that the model solution satisfies the model equations exactly.
This non-sequential method is derived from the sequential EnKF scheme by including the time dimension in the observation vector (Sakov et al., 2010). Consequently, all observations in the optimization period are gathered in a single observation vector y o , with an observation error covariance matrix R.
The optimization vector x, called the state vector in the classical method, contains all forcing and model parameters to be improved.
An ensemble of forcing x l , where l is the index of a member in the ensemble, is generated by adding a perturbation to the background (or original) forcing x b . If perturbed forcings are time dependent, their temporal evolution is gathered in a single vector, in the same way as observations.
For each perturbed forcing, the model is run. Note that at this step ensemble members are not influenced by observations. For each ensemble member, the observed part of the model state h(x l ) is extracted using a non-linear observation operator h(.). Consequently, model dynamic is included in the observation operator. Each element of the ensemble h(x l ) can be directly compared to the corresponding one in the observation vector y o .
Error covariance matrices can be estimated from the use of an ensemble of forcing state and model state. Using the same notation as in Barth et al. (2010Barth et al. ( , 2011, we define the deviation of the ensemble members around the ensemble mean of forcing S and observed part of surface currents E as where L is the number of ensemble members, the index l refers to the ensemble member (l = 1,. . . ,N) and . is the ensemble average operator. These matrices are scaled in such a way that SE T and EE T represent, respectively, ensemble forcing covariance and covariance between ensemble forcing and extracted observations: Note that there is no need to form explicitly and to store the covariance matrices, because analysis is computed in the ensemble space (e.g., Nerger et al., 2005). The optimal forcing x a that minimizes the analysis error is computed using Eq. (2), which is the Kalman filter analysis step with a non-linear observation operator (Chen and Snyder, 2007). It can be developed as or using the Sherman-Morrison-Woodbury formula: The superscripts "a" and "b" refer to the analysis and the background estimates. As in Kalman filtering, background and observation errors are assumed to be unbiased, but in the case of a strongly non-linear observation operator, such assumptions are not sufficient to ensure that the analysis x a is unbiased. However, the analysis is unbiased if Eq. (3) is verified (Barth et al., 2011). The model trajectory is corrected in a final step that consists in running the model with the optimized forcing x a . Consequently, observations influence the analyzed model solution using an optimal combination of ensemble forcing.