3.1 OMA and DINEOF

Some methods, i.e., open-boundary modal analysis (OMA) and data interpolating empirical orthogonal functions (DINEOFs), are nowadays widely used to fill spatiotemporal gaps in HFR measurements. We briefly explain these two methods below.

Open modal analysis (OMA, Kaplan and Lekien, 2007) is based on a set of linearly independent modes that are calculated before they are fit to the data. These modes describe all possible current patterns inside a two-dimensional domain (taking into account the open boundaries and the coastline). The amplitude of those modes is then fitted to current measurements inside the domain. OMA considers the kinematic constraints imposed on the velocity field by the coast since OMA modes are calculated taking into account the coastline by setting a zero normal flow. Depending on these constraints they can be limited in representing localized small-scale features as well as flow structures near open boundaries. Difficulties may arise when dealing with gappy data, especially when the horizontal gap size is larger than the minimal resolved length scale (Kaplan and Lekien, 2007) or when only data from one antenna are available. In the case of large gaps, unphysically fitted currents can be obtained if the size of the gap is larger than the smallest spatial scale of the modes, since the mode amplitudes are not sufficiently constrained by the data (Kaplan and Lekien, 2007). This is why when using OMA it is recommended to reach a compromise between the number of modes used for spatial scales larger than the largest gap and a sufficient number of modes correctly representing the spatial variance of the original fields (knowing that the spatial smoothing increases as the number of modes decreases). For this work, we used the OMA modules in the HFR Progs MATLAB package (https://cencalarchive.org/ProgsRealTime/HFR_Progs/, last access: 20 August 2018) to process radial velocities into total currents. Setting a minimum spatial scale of 20 km, 85 OMA modes were built for the regular grid shown in Fig. 1. OMA is especially useful to obtain a solution in the baseline area (i.e., the area between the two radar antennas for which total currents cannot be computed directly from radial information) and has been used to provide fields to compute accurate trajectories (Hernández-Carrasco et al., 2018; Solabarrieta et al., 2016) and surface Lagrangian transport in coastal basins (Hernández-Carrasco et al., 2018; Rubio et al., 2018).

The data interpolating EOF (DINEOF) is an iterative methodology used to interpolate gaps or missing data in geophysical datasets (Alvera-Azcárate et al., 2005; Beckers and Rixen, 2003). The technique is applied to an N×M data matrix, where N is the size of spatial locations of the geophysical field, and M the time dimension of the data. Before the methodology is initialized, part of the initially non-missing data is removed from the data matrix and stored for cross-validation. Then anomalies are computed by removing the time mean at each location, and gaps or missing data are replaced by zero-value anomalies. At this point the data matrix is iteratively decomposed and rebuilt by means of an EOF analysis with a fixed number of EOFs. Values of originally missing gaps evolve through this iteration until convergence is met. By changing the number of EOFs, a set of such iterations is produced using a different number of EOFs in each repetition. The optimal number of EOFs is then deduced by cross-validation comparing the initially stored data with the different versions of the reconstructed data. Once the optimal number is deduced, the originally retained data are introduced back in the data matrix, and a final reconstruction is made using the optimal number of EOFs.

DINEOF was introduced by Beckers and Rixen (2003) and Alvera-Azcárate et al. (2005, 2007) in both univariate and multivariate forms. The method was initially applied to sea surface temperature and ocean color (Alvera-Azcárate et al., 2005, 2015; Beckers et al., 2014; Esnaola et al., 2012; Ganzedo et al., 2011; Volpe et al., 2012), but the technique has also been used for other variables like sea surface salinity (Alvera-Azcárate et al., 2016). The DINEOF Fortran code (http://modb.oce.ulg.ac.be/mediawiki/index.php/DINEOF, last access: 20 August 2018) was applied to the combination of radial currents from the two antennas. The MATLAB–Octave complementary sources provided with the Fortran codes are also used to produce cross-validation masks internally required by DINEOF, following the procedure proposed in Alvera-Azcárate et al. (2005). As in the case of the above mentioned DINEOF applications, spatial locations with more than 95 % of missing data are not included in the DINEOF reconstruction procedure to prevent a negative impact on the quality of the overall reconstruction. The covariance filtering option proposed in Alvera-Azcárate et al. (2009) was also applied, setting the number of iterations of the filter to 12. The number of retained EOFs was high for all experiments (a minimum of 75) independently of the considered gap scenario (see Sect. 4).

3.2 SOM

Here a new methodology has been developed in order to reconstruct HF radar velocity fields in the statistical framework of the self-organizing map (SOM) analysis. SOM is a powerful visualization technique based on an unsupervised learning neural network, which is especially suited to extract patterns in large datasets (Kohonen, 1982, 1997). SOM is a nonlinear mapping implementation method that reduces the high-dimensional feature space of the input data to a lower-dimensional (usually two dimensional) network of units called neurons. In this way, SOM is able to compress the information contained in a large amount of data into a single set of maps. The learning process algorithm consists of a presentation of the input data to a preselected neuronal network, which is modified during an iterative process. Each neuron (or unit) is represented by a weight vector with the number of components equal to the dimension of the input sample data. In each iteration the neuron whose weight vector is closest (more similar) to the presented sample input data vector, called the best-matching unit (BMU), is updated together with its topological neighbors located at a distance less than the neighborhood radius R_n towards the input sample through a neighborhood function (see Hernández-Carrasco and Orfila, 2018, for a schematic representation of the SOM). Therefore, the resulting patterns will exhibit some similarity because the SOM process assumes that a single sample of data (input vector) contributes to the creation of more than one pattern, as the whole neighborhood around the best-matching pattern is also updated in each step of training. It also results in a more detailed assimilation of particular features appearing on neighboring patterns if the information from the original data enables it to do so. At the end of the training process, the probability density function of the input data is approximated by the SOM and each unit is associated with a reference pattern, with a number of components equal to the number of variables in the dataset, so this process can be interpreted as a local summary or generalization of similar observations.

For typical remote sensing imagery, SOM can be applied to both the space and time domain. Here, since we are interested in the reconstruction of HF currents, we have addressed the analysis in the spatial domain. In this case the input row vector has been built using the radial velocity maps at each time, so each neuron corresponds to a characteristic radial velocity spatial pattern over the coverage area of the HFR. Since each step iteration has an associated time and location of the sample, we can obtain the time of a particular spatial pattern computing the BMU for each time, providing a time series of the corresponding spatial pattern.

From these ideas, we can deduce the following simple algorithm for reconstructing missing values in the HFR velocity field from the available HFR data.

• i.

  Initialization (setup of the initial neural network). Each neuron (or SOM unit) has an associated weight vector, W, composed of random values of HFR radial velocities (called random initialization).

• ii.

  Training process. Radial HFR velocity values at a time t_i (Map of HFR velocities, U(x,t_i)) are used as input vectors, ${\mathbf{V}}_{\mathrm{i}}=\mathbf{U}(\mathit{x},t_i)$, to iteratively feed the neural network. This process is divided into two phases: a rough training using a larger neighborhood radius R_n in the neighboring function and a fine-tuning training for small R_n. Note that V_i values are the HFR velocity maps with missing points.

• iii.

  Spatial patterns. Once the training process has finished we obtain an ordered neural network of spatial patterns of radial HFR velocities (V_s). Now V_s does not include missing values.

• iv.

  Identification of missing values (NaNs). Locate the position and time of the HFR velocity grid points with missing values (x_m,t_m).

• v.

  Find BMU. Identify the most similar SOM spatial pattern to the HFR radial velocity map with missing values BMU(V_i(x_m,t_m)).

• vi.

  Replace missing values. Once we have assigned the BMU we replace the points with missing values with values of the corresponding grid point of the BMU, obtaining the filled HFR velocity field, V_f, where

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{V}}_{\mathrm{f}}(\mathit{x},t)=\mathrm{BMU}\left({\mathbf{V}}_{\mathrm{i}}\right({\mathit{x}}_{\mathrm{m}},t_{\mathrm{m}}\left)\right)\\ \text{(1)}& {\displaystyle }& {\displaystyle }\forall (\mathit{x},t),\mathbf{U}(\mathit{x},t)=\mathrm{NaN},\text{(2)}& {\displaystyle }& {\displaystyle }{\mathbf{V}}_{\mathrm{f}}(\mathit{x},t)={\mathbf{V}}_{\mathrm{i}}(\mathit{x},t)\forall (\mathit{x},t),\mathbf{U}(\mathit{x},t)\ne \mathrm{NaN}.\end{array}$$
• vii.

  Replace all NaNs. Repeat (v.) and (vi.) for all the points with missing values.

The ability of this method relies on the precision of identifying the proper BMU that accurately describes the missing dynamics. We try to optimize the algorithm by using as an input vector a concatenation of three maps of HFR velocities on three different dates. Thus we force the method to distinguish between three maps instead of one, avoiding the selection of a bad BMU, in particular when the HFR velocity map has a large number of missing points at time t₀ and there are not enough points to compare with the SOM patterns. The input vector is thus ${\mathbf{V}}_{\mathrm{i}}=\left[\mathbf{U}\right(\mathit{x},t_{\mathrm{0}}-\mathrm{\Delta }t),\mathbf{U}(\mathit{x},t_{\mathrm{0}}),\mathbf{U}(\mathit{x},t_{\mathrm{0}}+\mathrm{\Delta }t\left)\right]$. We have checked different values of Δt and found the best results when using Δt=4 h.

The initialization, training process and final output of the SOM algorithm have to be tuned in order to optimize the results and the computational cost by selecting particular control parameters. For instance, the optimal size of the neural network (number of neurons) depends on the number of samples and the complexity of the patterns to be analyzed. We choose the map size as [50×50], with 2500 neurons. Using different sizes, for instance [100×100], the spatial patterns are more detailed and more variability of patterns emerges. However, these new patterns make it difficult to find the proper BMUs for the specific date, creating considerable problems in accurately reconstructing the desired field in the missing points. The opposite happens using a reduced number of neurons: patterns are concentrated together in a few rough neurons without discriminating regions with different dynamical processes. We use a hexagonal map lattice to have equidistant neighbors and do not introduce anisotropy artifacts. Concerning the initialization, we opted for random mode. There are other initializations, i.e., based on a linear combination of the EOF modes of the HFR velocities, but all of them yield similar results. We choose random initialization since it is faster and missing data are accepted. For the training process we use the imputation batch training algorithm (Vatanen et al., 2015) adapted for data with missing values and a “Gaussian” type neighborhood function since this parameter configuration produces a good compromise between quantitative and topological error and computational cost (Liu et al., 2006). The neighborhood radius is chosen to vary from 20 to 0.1. Other values of the neighborhood radius have been tested without improving the results with the same computational cost. These SOM computations have been performed using the MATLAB toolbox of SOM v.2.0 Vesanto et al. (2000) provided by the Helsinki University of Technology (http://www.cis.hut.fi/somtoolbox/, last access: 20 August 2018).

Compared to conventional statistical methods like EOF and K-means, SOM is able to introduce nonlinear correlations and it does not require any particular functional relationship or distribution assumptions about the data, i.e., distribution normality or equality of the variance (Liu et al., 2006, 2016). SOM has been applied in a wide range of scientific disciplines. Among others, it has been used in climate sciences (Cavazos et al., 2002), genetics for DNA sequences (Nikkilä et al., 2002) and ecological sciences applications (Chon, 2011). In the physical and biological oceanography context SOM has also been used in several studies (Basterretxea et al., 2018; Charantonis et al., 2015; Hales et al., 2012; Hernández-Carrasco and Orfila, 2018; Liu and Weisberg, 2005; Liu et al., 2016; Richardson et al., 2003). However, to our knowledge, applications of SOM analysis to the reconstruction of HF radar velocity fields have not been addressed.

Figure 2Example of the four prototypes of KMA groups of gap distribution scenarios obtained from the KMA analysis applied to the HFR availability–absence matrix for 2014. Panel (a) represents a good data coverage scenario. Panel (b) accounts for the failure of one of the two antennas. Panel (c) characterizes a range decrease in both antennas. Panel (d) represents bearing coverage decreasing in one of the antennas. Matxitxako radials are plotted in red and Higer radials are plotted in blue (see Fig. S1 in the Supplement for the complete KMA results lattice).

Download

6.1 Eulerian comparison

Figure 5Scatterplot of the SOM (a), DINEOF (b) and OMA (c) reconstructed zonal velocities vs. observations for the four experiments (rows).

Download

Figure 6Scatterplot of the SOM model (a), DINEOF (b) and OMA (c) reconstructed meridional velocities vs. observations for the four experiments (rows).

Download

Scatterplots of zonal (u) and meridional (v) velocities reconstructed by the three models vs. the reference field and for the four experiments are shown in Figs. 5 and 6, respectively. In the case of OMA the cloud of points is, in general, wider with more points separated from the baseline. In general the behavior of the SOM and DINEOF methods is very similar in terms of RMS, BIAS and R² as displayed in Tables 1 and 2 for the u and v component of the velocity, respectively. Although SOM tends to underestimate currents and DINEOF to overestimate, OMA shows the worst reconstruction, especially for experiment C. It is worth noting that models behave similarly for component u and component v for all experiments. The lowest mean error is found for experiment D with values of 0.36, 0.35 and 0.38 for SOM, DINEOF and OMA, respectively, and the highest for experiment C with values of 0.43, 0.38 and 0.57 for SOM, DINEOF and OMA, with antenna failure as the worst-case gap scenario to be properly reconstructed. In particular, the scatterplot of reference vs. reconstructed zonal velocities from DINEOF in the case of the experiment C (Fig. 5 second column, third row) shows many points out of the main cloud. It suggests that the gap scenario C, in which a very small number of observations is available (failure of one antenna), is the most difficult gap scenario for the DINEOF methodology. Typically a 5 % minimum threshold of available data is used in most DINEOF applications reported in the literature to accept a candidate for a potential reconstruction (for instance, a gappy satellite image or an hourly HFR current field). As described in Sect. 3, such a threshold is also applied here as a preprocessing step before the technique is applied. The spatial field of the missing antenna is reconstructed using available data from the other antenna, as well as data from previous and subsequent time steps. Although in general the field might be acceptably reconstructed, it seems that some parts of the HFR image are not properly reconstructed. This result suggests that DINEOF should not be considered in situations represented in scenario C, especially when the antenna failure is persistent over time. Relative errors are slightly lower for DINEOF than for SOM but with more outliers as shown in Figs. 5 and 6.

Table 1Root mean square (RMS) error , mean bias (MB) and correlation coefficient (R²) of the zonal velocities (u) obtained using reconstructed HFR field from the three methodologies and for the four experiments with respect to the zonal velocities derived from the reference HFR currents.

Download Print Version | Download XLSX

Table 2Root mean square (RMS) error, mean bias (MB) and correlation coefficient (R²) of the meridional velocities (v) obtained using reconstructed HFR datasets from the three methodologies and for the four experiments with respect to the zonal velocities derived from the reference HFR currents.

Download Print Version | Download XLSX

6.2 Comparison of trajectories

Synthetic trajectories are computed by advecting particles in the original HFR velocity field (reference trajectories) and also in the OMA, SOM and DINEOF gap-filled currents. A total of 868 particles are initially uniformly distributed over the HF radar grid with an initial distance between particles of 5 km. Lagrangian particles are released every hour from 2 to 26 April 2012, 2013 and 2014 and advected during 72 h. Trajectories are computed using a fourth-order Runge–Kutta integration scheme and a bilinear interpolation of the gridded velocity field in space and linear in time. The mean distance of separation (D) between particle trajectories advected by the reference velocity field and by the reconstructed current fields averaged over all pairs of trajectories is plotted against time for each experiment in Fig. 7. After 10–15 h, differences in D start to be visible in all the cases. In all the experiments, DINEOF and SOM show lower separation distances and a similar behavior, while the OMA method shows the highest errors, in particular in the case of failure of one of the antennas (experiment C). In this case OMA capability to fill gaps is limited (as discussed in Kaplan and Lekien, 2007) and the effects can be noticed from the beginning of the trajectory simulation. The DINEOF and SOM methods show similar results, with separation distances around 1 km after 24 h of simulation and 3 km after 72 h. Values of D for SOM are slightly smaller in all of the experiments. In contrast, this distance reaches values of 5 km after 24 h in the case of OMA. The values of D obtained in our computations after 24 h of simulation are smaller than those reported in Molcard et al. (2009), Bellomo et al. (2015) and Kalampokis et al. (2016) for separation distances between real and virtual drifter trajectories. It is important to keep in mind that in these experiments the trajectories are calculated with data from different sources (HFR and drifters) and different samples (spatial resolution), which are collected in a different region (Mediterranean Sea). However, the fact that we obtain smaller values of D is an indication of the high performance of the three methods, with better results observed for the SOM method compared with OMA and DINEOF.

Figure 7Separation distance between particle trajectories advected by the filled HFR velocity field and by the reference velocity field as a function of time and averaged over all the pairs of trajectories. Different limits in the vertical axis of the figure have been used for a better data representation. The error bar is the confidence interval (1 standard deviation) of the spatial average of D(t).

Download

The spatial distribution of separation distances between real and simulated trajectories has also been analyzed in order to detect any anisotropy in the differences between them (Fig. 8). The results show the same behavior of the temporal separation distances between the different methodologies and experiments, but DINEOF and SOM yield better results than OMA. It is worth noting that OMA capability decreases in the western part of the study area in all the experiments, especially when there is a failure in one of the two antennas (experiment C). It is also noticeable that the time average of the separation distances is smaller in this experiment for DINEOF than for SOM. In general, all the methodologies show lower separation distances in the central and eastern part of the study area than in the western part, which demonstrates that the performance of the methodologies has spatial variations and is not uniform.

6.3 Lateral stirring analysis

Furthermore, we analyze the effect of the different gap-filling methodologies on the FSLE computations. FSLE was originally introduced in dynamical system theory to characterize the growth of non-infinitesimal perturbations in turbulence (Aurell et al., 1997) and to characterize flows with multiple spatiotemporal scales (Boffetta et al., 2000; Lacorata et al., 2001). In addition, FSLE can be used to unveil dynamical structures, called Lagrangian coherent structures (LCSs), thereby identifying fronts, eddies and jets (Hernández-Carrasco et al., 2011), as well as studying the flow stretching and contraction properties in geophysical data (Hernández-Carrasco et al., 2012; Joseph and Legras, 2002; Lipphardt et al., 2006; Shadden et al., 2009).

Figure 8Time average of the separation distance between trajectories computed from the filled HFR using the three methodologies and reference HFR velocities initiated in the same pixel.

Download

Figure 9Comparison of scale-dependent FSLE curves computed from the different filled HFR velocities averaged over 240 virtual particle-pair deployments homogeneously distributed through April 2012, 2013 and 2014. (a) Experiment A, (b) experiment C.

Download

The contribution of the scales captured by HFR in ocean surface dynamics can be analyzed by computing particle dispersion at different spatial scales, δ (Boffetta et al., 2000; Corrado et al., 2017; Haza et al., 2010), using the following equation for FSLE:

where 〈τ(δ,rδ)〉 is the time (τ) needed for the initial perturbation δ (separation of the initial pair of particles) to grow rδ averaged over all the pairs of particles for every initial perturbation δ and a fixed threshold rate, r. A small value, r<2, allows for the capturing of the relative dispersion driven by small coherent features, but r should not be close to 1 to avoid aliasing problems associated with the time step of particle advection at small scales (Haza et al., 2008; Lacorata et al., 2001). FSLE is a measure of relative dispersion for which the spatial variable is the independent variable.

We evaluate how the gap-filling methodologies impact the dynamical scales captured by the HFR. Figure 9 shows the averaged FSLE curves obtained from the three methodologies in two of the four experiments. The other experiments yield similar results. In all of them, λ(δ) is not constant over the range of scales covered by the HFR (scale dependent). In the REF case the best fit of the averaged FSLE curve shows a slope of −0.50 and a correlation coefficient, R², of 0.97. This indicates that the average scaling law of relative dispersion over the scales from 6 to 40 km is associated with a turbulent dispersion for flows with a k⁻² type of spectrum (Capet et al., 2008; Klein et al., 2008), suggesting that ageostrophic velocities and frontogenesis are contributing to lateral dispersion at these scales (Callies and Ferrari, 2013). This means that the separation rate is controlled by structures with sizes comparable to the separation itself, and therefore the dispersion regime is local. The slope of the FSLE curve derived from the REF HFR velocity field is in agreement with previous modeling studies, with a slope similar to those derived from surface synthetic drifter trajectories integrated by high-resolution numerical simulations (Choi et al., 2017). It further reinforces the fact that the velocity field of the HFR is adequate to study small-scale dynamics.

Comparing the FSLE slopes obtained from the three gap-filling methodologies with the REF HFR velocity field we find that the SOM methodology yields the most similar values, with slopes of $-\mathrm{0.54}\pm \mathrm{1}$ and R² greater than 0.96 for all the experiments, followed by the DINEOF with slopes of $-\mathrm{0.41}\pm \mathrm{1}$ and R² also greater than 0.96 for all the experiments. Finally OMA yields the largest differences with slopes between −0.27 and −0.30 and high R² (greater than 0.94). The SOM-FSLE slopes are slightly steeper than the REF-FSLE, suggesting that this methodology could introduce noise in the reconstruction of the data. On the other hand, DINEOF-FSLE slopes are slightly flatter than the REF-FSLE, which could be produced because this methodology smooths the velocity field or because it favors larger scales as captured by EOFs. However, OMA-FSLE slopes are clearly flatter than the REF-FSLE. In this case, λ(δ) is almost constant over the range of scales (scale independent), the separation of particles is almost exponential with a constant rate and we have a nearly exponential regime. In this regime the relative dispersion is nonlocal, indicating that not only small-scale processes govern the dispersion but also larger structures. This result suggests that OMA has a strong smoothing character, removing small features from the velocity field, as reported in previous Eulerian studies (Kaplan and Lekien, 2007).

Compared with other observational studies using HFR the λ obtained from our computations is 1 order of magnitude smaller than the one reported in Haza et al. (2010) (λ=4 to 7 days⁻¹ for δ<1 km) computed from VHF radar with 300 m of spatial resolution in the Gulf of La Spezia. This is expected when comparing their scales of interest (0.1–1 km) against our range of resolved scales (5–25 km) (Corrado et al., 2017). Hernández-Carrasco et al. (2018), using HFR velocities at 3 km of resolution in the Ibiza Channel, reported larger values of λ (between 0.6 and 0.2 days⁻¹) at spatial scales between 3 and 30 km following a Richardson slope. Other studies based on drifters in the Gulf Stream (Lumpkin and Elipot, 2010), the Gulf of Mexico (Poje et al., 2014) and the Mediterranean Sea (Griffa et al., 2013; Schroeder et al., 2011) have found higher values of λ(δ) at these scales. Regional differences and possible observational biases could explain these discrepancies in FSLE values. For instance, differences in the topographic boundaries and/or the presence of coastal jets, tidal currents or particular wind regimes can influence the dynamics affecting transport processes at the coastal sea surface.

In general the λ values obtained for the three methodologies and for all the experiments are smaller than for REF, in particular when OMA velocities are used. At the smallest scales studied here (5–15 km) SOM-FSLE is the most similar, and at larger scales DINEOF-FSLE and OMA-FSLE values are closer to REF-FSLE.

Next we analyze the effect of the reconstructed velocities on the LCS. FSLE is used to obtain the LCS by computing the minimum growth time of particle-pair separations from δ to rδ among the four neighbors for each position in order to obtain an FSLE map. In this case, r has to take large values (r≫2) to adequately distinguish regions of extrema in the FSLE field. Note that in this case the average over the pair of particles in Eq. (6) is omitted.

First we compare some snapshots of the LCS derived from the three reconstructed HFRs in order to see differences in the dynamical structures (Fig. 10). We only show LCS from two experiments, A and C, since the results obtained in experiments B and D are similar to those of experiment A. The computed Lagrangian flow structures look rather the same in the case of SOM and DINEOF, despite the large number of missing points reconstructed in both experiments. This is also confirmed by computing the histograms of the FSLE (not shown), which turn out very similar. This is likely due to the averaging process along trajectories performed to compute the FSLE, which tends to remove possible errors introduced in the reconstructed velocity fields (Hernández-Carrasco et al., 2011). FSLE is robust against the gaps in the HFR velocity field. On the other hand, OMA is the most affected methodology regarding LCS computations. Even if not very different, one can see some discrepancies in the shape and the location of OMA-LCS with respect to the REF-LCS, as happens in experiment C (Fig. 10b).

The time evolution of the spatial average over the HFR domain of the absolute relative error (ARE; Sect. 5) of the attracting LCS computed from the reconstructed velocity fields with respect to the REF-LCS is plotted in Fig. 11. A common characteristic for the three methods is large values of ARE during early April 2012 and middle April 2014 for all the experiments, in agreement with the time periods with the largest number of missing values. Comparing the methods one can see that the relative error is mostly greater for the case of OMA and for all the experiments, in particular during the second half of April 2012 and the middle of April 2013 and 2014. All the methods are more affected by the antenna failure simulated in experiment C, in accordance with the Eulerian comparison.

Figure 10Snapshots of attracting LCS computed from the three filled and the reference HFR currents for experiment A (a) and experiment C (b) corresponding to 11 April 2014, 00:00 UTC.

Download

.

Figure 11Time series of the spatial average of the pixel-by-pixel absolute relative error of the LCS computed from the three filled HFR velocities with respect to LCS obtained from the reference velocity field.

Download

Further analysis is conducted by performing a regional characterization of the impact of the gap-filling methods on the LCS looking at the spatial distribution of the time-averaged values of ARE (Fig. 12). One can appreciate that points with accumulated large errors are located in the same regions for SOM and DINEOF, i.e., in the southern area of the BoB HFR domain, and that these regions are the same for the four experiments (near the southern coastline). In addition to possible smoothing or noise, a discrepancy in the location of LCS with respect to the REF-LCS could explain the large values of the errors. In the case of OMA, additional regions with large errors are distributed all over the BoB HFR domain.

Finally, we quantify the difference between the LCS computed from the filled velocities and the REF-LCS by computing the spatial and temporal mean of the relative error defined in Eqs. (3) and (4) over all the LCS maps. We summarize in Table 3 the results of these computations. The most significant result is that ARE values are similar for SOM and DINEOF and for all the experiments (around 0.21), being in general slightly smaller for the case of SOM. The largest ARE is found for OMA-LCS in experiment C with the worst-case gap scenario affecting LCS, with a value of 0.27. The values of MB show that in general the LCSs are slightly overestimated. Only in experiment C are the values of OMA-LCS underestimated.

Figure 12Maps of the time average of the pixel-by-pixel absolute relative error of the LCS computed from the three filled velocity field with respect to LCS obtained from the reference currents. White denotes points at which LCSs have not been identified.

Download

Our results show that even when a large number of pixels is missing the FSLEs still give an accurate picture of oceanic transport properties, and the velocity field filled by the three methods analyzed does not introduce artifacts in the Lagrangian computations.

6.4 Residence times

Another Lagrangian quantity suitable to describe transport process is residence time (RT) (Buffoni et al., 1996; Hernández-Carrasco et al., 2013; Lipphardt et al., 2006). RT is commonly used to characterize the fluid interchange between different oceanic regions and is defined as the interval of time that a particle remains in an area before crossing a particular boundary. In the calculation of RT particles are initially located on the grid points of the HFR velocity field and are integrated in time during 14 days. We consider only 14 days owing to the short period of available data without missing values (REF field). In these computations we assign the maximum possible value of RT (14 days) to the particles that remain in the area without crossing the preselected boundary after the 14 days of integration.

Table 3Mean absolute relative error (ARE) and mean bias (MB) of the LCS obtained using reconstructed HFR velocities from the three methodologies with respect to the REF-LCS for the four experiments.

Download Print Version | Download XLSX

As done for the LCS, a first comparison is performed looking at the spatial distribution of RT obtained from the three different filled HFR datasets. In general, RT maps from the three methodologies are similar to the REF-RT (maps of RT obtained from the reference HFR velocity field). As an example, Fig. 13 shows two maps of RT for experiments A and C corresponding to 4 April 2013, 18:00 UTC. In this figure we color the initial position of particles in the BoB according to the time they need to escape from the HFR domain. Regions with particles with short residence times are shown in cyan–blue, and initial positions of particles having longer residence times are indicated in light brown. Again these two experiments are representative for the other experiments since experiments B and D produce an RT spatial distribution similar to that shown for experiment A (Fig. 13a). As it can be appreciated, the spatial patterns of RT obtained from the filled HF radar velocity fields are similar to the REF-RT in experiment A (and also for experiments B and D, not shown). Only some differences are found in OMA-RT maps compared to REF-RT for experiment C, as can be seen in Fig. 13b (OMA panel) with higher values of RT in the west and southeast regions of the HFR domain.

To further analyze the periods when the gap-filling methods introduce more errors in the RT computations we plot in Fig. 14 the time series of the spatial averages of the ARE of the RT. Most of the reconstructed derived RT values approximate REF-RT values, with ARE smaller than 0.2 along the time series. It can be seen that, except at the beginning of April 2012, values of ARE for OMA-RT are larger than SOM-RT and DINEOF-RT, in particular at the end of April 2012 and in the middle of April 2013 and April 2014. Again, experiment C turns out as the more critical situation for gap distribution to be reconstructed, reaching the largest values of ARE (around 0.25).

Figure 13Snapshots of RT computed from the three filled velocity fields and the reference HFR currents for experiments A and C corresponding to 4 April 2013, 18:00 UTC.

Download

Figure 14Time series of the spatial average of the pixel-by-pixel absolute relative error of the RT computed with the reconstructed velocities with respect to the REF-RT.

Download

Figure 15Maps of the time average of the pixel-by-pixel absolute relative error of the RT computed from the three filled velocities with respect to RT obtained from the reference velocity field.

Download

In order to reveal regions where the effect of the reconstruction on the RT is higher we compute the time-averaged relative error for all experiments (Fig. 15). This figure shows that, as for the LCS comparison, regions of accumulated errors are the same for SOM and DINEOF and this spatial distribution is quite similar for the four experiments. OMA shows different regions with large errors, being the method with the greatest number of pixels with large error. A common feature is that the three methods yield large errors near the southern coastline and in the west and northwest area of the HFR domain. It suggests that the three methods remove some dynamical features in this region. For points located in central regions (around 2^∘ W and 44^∘ N) values of RT derived from the reconstructed velocities better approximate the REF-RT values. The greater availability of data over this region could explain this agreement in the RT computations.

To finish this section, we summarize in Table 4 the results of the ARE and ME of RT averaged over all the RT values. Although ARE is slightly greater in the case of OMA than using SOM and DINEOF, in general these errors are low in all cases, with 0.16 as the maximum ARE obtained over all the experiments. The absolute error is smaller for SOM and DINEOF (between 0.12 and 0.13 for both cases) and slightly larger for OMA (between 0.13 and 0.16). Furthermore, the AREs for SOM-RT and DINEOF-RT are quite identical. In all cases the MBs are positive, meaning that RTs are overestimated. The largest values of MB are found, on average, for OMA-RT (MB ranging from 0.7 h for experiment D to 8.1 h for experiment A), followed by SOM-RT (MB ranging from 5.0 h for experiment B to 7.8 h for experiment C) and finally by DINEOF-RT (MB ranging from 4.3 h for experiment B to 5.4 h for experiment A).