Impact of HF radar currents gap-filling methodologies on the Lagrangian assessment of coastal dynamics

In coastal basins High-frequency radar, HFR, is a cost-effective monitoring technique that allows to obtain high-resolution continuous surface currents, providing new insights for understanding small-scale coastal ocean transport and dispersion processes. In the last years the use of Lagrangian diagnosis to study mixing and transport properties is growing in importance. A common condition among all the Lagrangian techniques is that complete spatial and temporal velocity data is required to 5 compute the trajectories of virtual particles in the flow. However, hardware or software failures in HFR can compromise the availability of data, resulting in incomplete spatial coverage fields or non-data periods. In this regard, several methods have been widely used to fill spatio-temporal gaps in HFR measurements. Despite the growing relevance of these systems there are still many open questions concerning the reliability of the gap-filling methods for the Lagrangian assessment of the coastal ocean dynamics. In this paper, we first develop a new methodology to fill gaps in the HFR velocity field based on Self-Organizing 10 Maps (SOM). Then a comparative analysis of this method with the most extended available methods for gap-filling in HFR systems, i.e., Open-Boundary Modal Analysis (OMA) and Data Interpolating Empirical Orthogonal Functions (DINEOF), is performed. The error introduced by each approach is quantified in the Lagrangian computations, i.e., finite size Lyapunov exponents, Lagrangian Coherent Structures and Residence Times. We determine the limit of applicability of each method regarding four experiments based on the typical temporal and spatial gaps distribution observed in HFR systems unveiled by a k-means 15 clustering analysis. The data sets used for this study are the hourly radial velocities obtained from the HFR system located in the SE Bay of Biscay, with a spatial resolution of 5 km in range and 5◦ in angle. Our, results, show that even when the number of missing data is large, the Lagrangian diagnosis still give an accurate description of the oceanic transport properties.


Introduction
Knowledge of the spatial and temporal complexity of the coastal processes is one of the major challenges in oceanography.Understanding the coastal hydrodynamics is crucial to quantifying the contribution of the coastal margins to the ocean's biological productivity (Pauly et al., 2008;Cloern et al., 2014) and to the global CO 2 storage (Bauer et al., 2013).
Large-scale currents in coastal regions are strongly affected by the intricate orography and by the influence of the local boundary layers: the input of energy at the surface through heterogeneous heat fluxes and local wind forcing and the dissipation at the coast and at the seabed.This results in a highly variable nearshore circulation characterized by a complex interaction of multiscale processes that requires specific numerical and observational approaches to properly study their variability (Haidvogel et al., 2000).High Frequency Radars (HFRs) data can potentially be used to observe and provide information of these processes through hourly measurements of the surface currents at high-resolution.
HFR systems were originally developed in Stewart and Joy (1974) and Barrick et al. (1977) to remotely sensing the state of the sea surface based on the back-scatter HF radio waves to the measurements of the ocean surface currents.
Over the past 20-25 years there has been a rapid growth in the use of coastal radars demonstrating the possibility of observing and monitoring the complex surface current dynamics, and leading to a rapid spread of installations of integrated HFR observatories in many coastal regions.Currently, HFR is the unique technology able to provide autonomously continuous hourly surface velocity measurements over wide coastal areas (typical range of 30-100 km from the coast) and at high spatial (a few km) resolution depending on the working frequency (Paduan and Washburn, 2013;Bellomo et al., 2015;Lana et al., 2016;Rubio et al., 2017).Contrary to satellite altimetry, the use of HFR data can potentially provide gridded velocities with the necessary resolution in both space and time to unravel the role of small scales on the dynamical properties (Hernández-Carrasco et al., 2018).
In particular, HFR observations are crucial in coastal areas for those applications associated with transport processes, not only addressing the study of marine ecosystems but also for a wide range of coastal activities.These include search and rescue operations (Ullman et al., 2006); predicting and mitigating the spread of oil spill or other pollutants (Lekien et al., 2005;Shadden et al., 2009); understanding the impact of transport and mixing properties on relevant biogeochemical properties such as the primary productivity of plankton (Cianelli et al., 2017;Hernández-Carrasco et al., 2018); coastal larval transport or fishery management (Bjorkstedt and Roughgarden, 1997).
Lagrangian approach addresses the effects of the velocity field on transported substances, which clearly is of utmost relevance for studying transport processes.These approaches, have the advantage of exploiting both spatial and temporal variability of a given velocity field.They can even unveil sub-grid filaments generated by chaotic stirring, providing a more complete description of transport phenomena.Specially, the concept of Lagrangian Coherent Structure (see the review by Haller (2015)) has been receiving growing attention in the last decades in geophysical sciences since it gives an unprecedented description of transport processes in a specific area, identifying local maximizing lines of material stretching, separating regions with different dynamical behavior, and detecting avenues and barriers to transport, which are of great relevance for the marine biological dynamics (Tew Kai et al., 2009;Hernández-Carrasco et al., 2014;Bettencourt et al., 2015).Lagrangian diagnostics, including LCS, residence time, Lagrangian divergence, provide an improved description of the flow geometry and variability, which enables the study of the effect of transport properties derived from HFR on biological quantities in coastal seas (Rubio et al., 2018;Hernández-Carrasco et al., 2018), to reduce the damaging effects of coastal pollution Lekien et al. (2005), as well as to localize zones which are relevant for the tracking of debris accumulation or jellyfish aggregations.
Lagrangian diagnosis require the complete spatial and temporal velocity data in a specific region to compute trajectories of synthetic particles.However, HFR fails either by external circumstances, i.e., environmental factors, or due to hardware or software malfunction, power failure and communication disruptions at the radar stations, leading to incomplete spatial measurements in form of gaps of data both in space and time.In general the occurrence of gaps in the data follows a pattern which can be associated with a particular cause, for instance the instrumentation malfunction, sea state, signal interference, antenna configuration, facilitating their interpretation and simulation.
In this regard, different methodologies have been developed in the last years with the aim of filling velocity fields derived from HFR measurements.The most widely extended technique is the Open-Boundary Modal Analysis (OMA) (Kaplan and Lekien, 2007).However other methods have been proved to be efficient, namely; DINEOF method based on empirical orthogonal functions (Alvera-Azcárate et al., 2005); including a functional penalizing the OMA modes variability (Yaremchuk and Sentchev, 2009); DCT-PLS method based on penalized least squares regression method using a three-dimensional discrete cosine transform (Fredj et al., 2016), among others.Although the performance of some of these methods has been analyzed from the Eulerian point of view, that is, comparing the original and the reconstructed velocity fields (Fredj et al., 2016;Yaremchuk and Sentchev, 2009;Kaplan and Lekien, 2007), an assessment of the reliability of the Lagrangian metrics when computed with reconstructed data is still lacking.
Here we introduce a new method based on unsupervised neural networks, Self-Organizing Maps (SOM).Then, we perform an intercomparison of this new methodology with OMA and DINEOF.We analyze the effect of these gap-filling methods on the Lagrangian diagnosis derived from HFR velocities in the SE of the Bay of Biscay (BoB).Their robustness under different scenarios of data gaps are discussed.
The paper is organized as follows.After this introduction, we first describe the data.Gap-filling methods used in the comparison, including a detailed description of the new method developed here are described in section 3. Section 4 is devoted to describe the different data gaps scenarios commonly found in the SE BoB HFR system and used in the evaluation.The results of the Eulerian and Lagrangian comparison are shown in Sec. 5. Department.The sites are located in Cape Higer and Cape Matxitako (Figure 1).The averaged Doppler backscatter spectrum obtained from the received signal (in a window of 3 h) is processed to obtain hourly radial currents using the MUSIC algorithm (Schmidt, 1986).The coverage of radial data is up to 150 km with a 5 km range cell resolution and 5º angular resolution.Radial data are quality controlled using advanced procedures based on velocity and variance thresholds, noise to signal ratios and radial total coverage.Successful validation analysis for this HFR with surface drifters and mooring data have been previously done in L. et al. (2014) and Rubio et al. (2011).Hourly HFR data for April 2012, April 2013 and April 2014 have been used in the HFR system of the SE BoB. 3 Gap-filling methods 3.1 OMA and DINEOF Some methods, i.e., Open-boundary Modal Analysis (OMA); Data Interpolating Empirical Orthogonal Functions (DINEOF), are nowadays widely used to fill spatio-temporal gaps in HFR measurements.We briefly explain these two methods.

HF radar data
The 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 and 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 fit to current measurements inside the domain.The OMA analysis considers the kinematic constraints imposed on the velocity field by the coast since the OMA modes are calculated taking into account the coastline by setting a zero normal flow constraint.Depending on these constraints they can be limited in representing localized small-scale features as well as flow structures near open boundaries.Besides, 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 there is just data from one antenna.In the case of large gaps, unphysically fitted currents can occur in areas without data if the gaps are larger than the smallest spatial scale of the modes, since the mode amplitudes are not sufficiently constrained by the data (Kaplan and Lekien, 2007).That is why for the practical application of OMA it is recommended to get a compromise between limiting the number of modes used to those whose spatial scale is larger than the largest gap and allowing a sufficient number of modes to correctly represent 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/ocmpmb/COCMPwiki/) to process radial velocities into total currents.Setting a minimum spatial scale of 20km, 85 OMA modes were built for the regular grid shown in Figure 1.OMA is especially useful to obtain a solution in the base line area (i.e. the area between the two radar antennas where total currents cannot be computed directly from radial information) and has been used to provide fields to compute accurate trajectories (Solabarrieta et al., 2016;Hernández-Carrasco et al., 2018) and surface Lagrangian transport in coastal basins (Rubio et al., 2018;Hernández-Carrasco et al., 2018).
The Data Interpolating EOFs (DINEOF) is an EOF-based iterative methodology used to interpolate gaps or missing data in geophysical datasets (Beckers and Rixen, 2003;Alvera-Azcárate et al., 2005).The technique is applied to a 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 a convergence is met.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.
The DINEOF Fortran code (http://modb.oce.ulg.ac.be/mediawiki/index.php/DINEOF) 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 by Alvera-Azcárate et al. (2005).Like in the case of the above mentioned DINEOF applications, spatial locations with more than a 95% of missing data are not included in the DINEOF reconstruction procedure to prevent to negatively impact the quality of the overall reconstruction.The covariance filtering option proposed in Alvera-Azcárate et al. ( 2009) was also applied using the standard filter options proposed by these authors and 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 Sec. 4).

SOM
Here a new methodology has been developed in order to reconstruct HF radar velocity fields in the statistical framework of the Self-Organizing Maps (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(Kohonen, , 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) networks of units called neurons.In this way, SOM is able to compress the information contained in a large amount of data in a single set of maps.The learning process algorithm consists in 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 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.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 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, the SOM can be applied to both space and time domain.Here, since we are interested in the reconstruction of HF currents, we have addressed the analysis in the spatial domain of the datasets.In this case the input row vector has been built using the radial velocities maps at each time, so each neuron corresponds to a characteristic radial velocities spatial pattern, composed of a specific spatial distribution of different values of radial velocities over the coverage area of the HFR.Since each step iteration has associated a 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: Set up of the initial neural network.Each neuron (or SOM unit) has associated a weight vector, W, composed of random values of HFR radial velocities (called Random initialization).
(ii) Training process: Radial HFR velocities values at a time t i (Map of HFR velocities, U(x,t i )) are used as input vectors, V i =U(x,t i ), to iteratively feed the neural network.This process is divided in two phases: a rough training using a larger neighborhood radius R n in the neighboring function and a finetune training for small R n .Note that V i are the HFR velocity maps with missing points.
(iii) Spatial patterns: After the training process 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 vales (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 (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 describes accurately the missing dynamics.We try to optimize the algorithm by using as input vector a concatenation of three maps of HFR velocities at three different dates.Thus we force the method to distinguish between three maps instead of one, avoiding the selection of a bad BMUs, in particular when the HFR velocity map has a large number of missing points at the time t 0 and there are not the enough number of points to compare with the SOM patterns.The input vector is thus have checked different values of ∆t and we have found that the best results are found when using ∆t = 4 hours.
Initialization, training process and final output of the SOM algorithm have to be tuned up in order to optimize the results and the computational cost by selecting particular control parameters.For instance, the optimal size of map (number of neurons) depends on the number of samples and on the complexity of the patterns to be analyzed.We choose the map size as [50 x 50], with 2500 neurons.Using different sizes, for instance [100 x 100], the spatial patterns are more detailed and more variability of patterns emerge.However, these new patterns make difficult to find the proper BMUs for the specific date providing considerable problems to accurately reconstruct the desire field in the mixing points.The opposite happens using a reduced number of neurons, patterns are concentrated together in few rough patterns without discriminating regions with different dynamical processes.We use hexagonal map lattice in order to have equidistant neighbors and do not introduce anisotropy artifact.Concerning to the initialization, we opted for random mode.There are other initializations, i.e. based on linear combination of the EOFs modes of the HFR velocities, but all of them yield similar results.We choose random initialization since it is faster as well as missing data is accepted.For the training process we use the imputation batch training algorithm (Vatanen et al., 2015) adapted for data with missing values, and '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/).
Comparing with the conventional statistical methods like the 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(Liu et al., , 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), in genetics, working with DNA sequences (Nikkilä et al., 2002), or ecological sciences applications (Chon, 2011).In the physical and biological oceanography context SOM has also been used in several studies (Charantonis et al., 2015;Liu and Weisberg, 2005;Liu et al., 2016;Richardson et al., 2003;Hales et al., 2012;Hernández-Carrasco and Orfila, 2018b).However, to our knowledge, applications of SOM analysis on the reconstruction of HF radar velocity fields have not been addressed.

Experiments: Catalog of HFR Gaps
The most common real scenarios for spatial gaps in HFR data are mainly represented by individual antenna failures, range and/or bearing reduction.The radio signal emitted by an HFR travels along and back through the ocean surface due to the conductivity of the ocean and the currents velocity is measured based on the Bragg scattering phenomena (Barrick et al., 1977) of the received signal.Any affection to this process will result in gaps in the final data.Among the most common reasons to find spatio-temporal data gaps in radial (and total) are adverse environmental conditions and/or electromagnetic problems as the lack of Bragg scattering ocean waves, severe ocean wave conditions, the occurrence of radio interference or changes of the electromagnetic field in the vicinity of the antennas which lead to invalid antenna patterns and calibration parameters.
Additionally there is a permanent region between the antennas, the so-called baseline, where the total currents cannot be computed in an accurate way.The baseline between two HFR sites is defined as the area where the radial components from the two sites make an angle of less than 30°, so the total velocity vectors created from radial data within this data contain greater uncertainties.The solution in the baseline is normally not computed, so we observe a data gap in an area which is delimited by the rule of GDOP (Barrick, 2002;Chapman et al.) and the limits set to this quality control parameter in the processing of the data from radials to totals.The permanent spatial gap in the baseline, frequently located near the coastal area between the antennas, is an issue than can be problematic for the assessment of the dynamics of some areas.
In order to characterize the most typical and realistic gap types observed in the Basque HFR system, K-means classification algorithm (Hastie et al., 2009), henceforth KMA, is applied to the real radial data for 2014.Previous to applying the KMA, the radial data are converted to 1 or 0 values, depending on the availability or absence of data for each radial position, respectively.
Then KMA are used to classify the dataset into a specified number of groups according to the similarity in the distribution of gaps exhibited in the HFR data set (Hastie et al., 2009).The selection of the number of groups was done qualitatively, as proposed by Guanche (2014).In our case we choose to keep 16 groups since this number of groups facilitates the interpretation of the gaps scenarios without loosing variability of the data set (Figure S1   associated to the patterns with a percentage of gaps larger than a 30%, for this system, is likely due to a failure of one of the two antennas than to other causes.This last scenario represents a situation where data in selected locations have been removed (for instance after a quality control procedure based on a velocity or variance threshold, as recommended by the QARTOD manual, IOOS 2016).Figures 3 and 4 show the temporal and spatial distribution of the coverage % of the original data and the new dataset generated by introducing the previously described gaps to the original data.
(D) Random gaps: spatial gaps in 30introduced in time.K-means analysis shows that the situation associated to the patterns with a percentage of gaps larger than a 30to a failure of the antenna than to other causes.This last scenario represents a situation where data in selected locations have been removed (for instance after a quality control procedure based on a velocity or variance threshold, as recommended by the QARTOD manual, IOOS 2016).
These gap scenarios (from now on referred as experiment A, B, C and D) are used to test the SOM, DINEOF and OMA gap-filling methodologies evaluating their robustness from the Lagrangian point of view.For OMA, total maps are generated using OMA directly on radial data.For DINEOF and SOM, hourly radials are gap-filled first and then totals are generated using the same least mean square algorithm (spatial interpolation radius of 10km) used to build the reference data series (i.e. from the reference radial files with no gaps).

Definitions and statistical metrics used for the filled-original currents comparison
In this study we compute different Lagrangian quantities using the HFR velocity field filled by the three different methods We use conventional statistical metrics to measure difference between gap-filled and observation data to quantify methodology skill with respect to observed HFR data: absolute relative error (ARE), mean bias (MB), and root mean squared error (RMS).Denoting a set of reference observational values as R, corresponding data obtained from the filled velocity field as G, using brackets to denote the mean of the set, and a prime to denote perturbations from the mean, i.e., G =G-< G >, these error metrics are defined as: where N is the number of pixels in the velocity field.

Eulerian Comparison
Scatter plots of zonal (u) and meridional (v) velocities reconstructed by the three models against the reference field and for the , where a very small number of observations is available (failure of one antenna), is the most aggressive situation for 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, an hourly HFR current field, etc.).As described in Sec 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 being reconstructed using the available data from the other antenna, as well as with data from previous and posterior time steps.Although part of the field might be acceptably reconstructed, it seems that there are some parts of the HFR image that are not being properly reconstructed.This result suggests that DINEOF should not be considered in situations represented scenario C, specially when the antenna failure is persistent in time.Relative errors are slightly 10 lower for DINEOF than for SOM but with more outliers as shown in Figures 5 and 6.

Comparison of trajectories
Synthetic trajectories are computed advecting particles in the HFR original dataset (to be used as reference) and also in the OMA, SOM and DINEOF gap-filled currents to make comparisons.Trajectories are computed using a fouth-order Runge-Kutta integration scheme and a bilinear interpolation of the gridded velocity field both in space and time.The mean distance of separation averaged over all the comparison pairs of trajectories (reference trajectories vs. simulated trajectories with the 5 reconstructed current fields) has been plotted against time for each experiment in Figure 7.After 10-15 hours, the differences start to be visible in all the cases.In all the experiments, DINEOF and SOM show lower separation distances and a similar behavior, while OMA method is the one showing the highest errors, in particular in the case of failure of one of the antennas (Experiment C).In this case OMA capability to fill the gaps is limited (as discussed in Kaplan and Lekien (2007)) and the  Sea).However, the fact that we obtain smaller values of D is an indication that the performance of the three methods used here is very high, with better results observed for the SOM method as compared with OMA and DINEOF.
The spatial distribution of longitude and latitude distances between real and simulated trajectories has also been analyzed in 10 order to detect any anisotropy in the differences between them (Figure 8).The results show the same behavior of the temporal separation distances between the different methodologies and experiments; DINEOF and SOM yield better results than OMA.
It is worth to note that OMA capability decreases in the half west part of the study area in all the experiments and specially when there is a functioning failure in one of the two antennas (Experiment C).It is also noticeable that the separation distance values are lower in this experiment for DINEOF than for SOM.This was not appreciable in the time comparison of the distance of separation figures but we see it in the spatial distribution after 72 hours of integration.Generally, all the methodologies show lower separation distances in the central and East part of the study area than in the West which demonstrates that the capabilities of the methodologies have spatial variations and they are not uniform.

Lateral stirring analysis
Furthermore, we analyze the effect of the different gap-filling methodologies on the FSLEs computations.FSLE was originally introduced in the dynamical system theory to characterize the growth of non-infinitesimal perturbations in turbulence (Aurell et al., 1997) and to characterize flows with multiple spatio-temporal scales (Boffetta et al., 2000;Lacorata et al., 2001).In addition, FSLE can be used to unveil dynamical structures, called Lagrangian Coherent Structures, LCS, identifying fronts, eddies, jets, (Hernández-Carrasco et al., 2011), as well as to study the flow stretching and contraction properties in geophysical data (Joseph and Legras, 2002;Lipphardt et al., 2006;Hernández-Carrasco et al., 2012;Shadden et al., 2009).
The contribution of the scales captured by the HFR in the ocean surface dynamics can be analyzed by computing particle dispersion at different spatial scales, δ (Boffetta et al., 2000;Haza et al., 2010;Corrado et al., 2017)  equation for FSLE, where τ (δ, rδ) is the time (τ ) needed for the initial perturbation δ (pair-particles separation) to grow rδ averaged over all the particles pairs for every initial perturbation δ and a fixed threshold rate, r.A small value (r < 2) allows capturing the relative dispersion driven by small coherent features and not close to 1 to avoid aliasing problems associated with the time step of 5 particle advection at small scales (Lacorata et al., 2001;Haza et al., 2008).FSLE is a measure of relative dispersion where the spatial variable is the independent variable.We evaluate how the fact of filling the gaps affects the dynamical scales resolved 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 bestfitting of the averaged FSLE curve shows a slope of -0.50 with a correlation coefficient, R 2 , of 0.97.This indicates that the average scaling law of the relative dispersion over the scales from 6km to 40km is associated to a turbulent dispersion for flows with a k −2 type of spectrum (Klein et al., 2008;Capet et al., 2008), suggesting that ageostrophic velocities and frontogenesis are contributing in the lateral dispersion at these scales (Callies and Ferrari, 2013).This means that the separation rate is controlled by structures with size comparable with the separation itself, 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 slopes obtained from the three methodologies with the REF HFR velocity field we find that SOM methodology yields the most similar, with slopes of -0.54 ± 1, with a R 2 larger 0.96 for all the experimets, followed by the DINEOF with slopes of -0.41 ± 1, with R 2 also larger than 0.96 for all the experiments.Finally OMA yield the larger differences with slopes comprised between -0.27 and -0.30 with high R 2 (larger than 0.94).The SOM-FSLE slopes are slightly steeper than the REF-FSLE suggesting that this methodology might be introducing 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 is smoothing the velocity field or because it favors larger scales as captured by EOFs.However OMA-FSLE slopes are much flatter than the REF-FSLE.In this case, λ(δ) is almost constant over the range of scales (scale independent), the separation of particles is exponential with a constant rate and we have an exponential regime.In this regime the relative dispersion is non-local, that is not only the small scale processes govern the dispersion but also larger structures.It 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).2018), using HFR velocities at 3km of resolution in the Ibiza Channel, reported larger values of λ (between 0.6 and 0.2 days −1 ) 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) or in the Gulf of Mexico (Poje et al., 2014) or in the Mediterranean Sea (Schroeder et al., 2011;Griffa et al., 2013) have found higher values of λ(δ) at these scales.Regional differences as well as possible observational biases, could explain these discrepancies values of FSLE.
For instance, differences in the topographic boundaries, the presence of coastal jets, tidal currents, particular wind regimes, etc, can influence the dynamics affecting the transport processes at the coastal sea surface.
In general the λ 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 evaluate the effect on the LCS obtained from the three methodologies with the REF-LCS.FSLE is used to obtain the LCS by computing the minimum growth time of pair-particles separations from δ to rδ among the four neighbors for each position in order to obtain a 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 HFR in order to see how different can be the dynamical structures (Fig. 10).We only show LCS from two experiments, A and C, since the results obtained in the experiments B, and D are similar to those of the experiment A. The computed Lagrangian 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 see some differences in the shape as well as in the location of OMA-LCS with respect to the REF-LCS, as happens in the Experiment C (Fig. 10 b).Further analysis is addressed performing a regional characterization of the impact of the gap-filling methods on the LCS looking at the spatial distribution of the time average values of ARE (Fig. 12).One can appreciate that points with accumulated large errors are located in the same regions for SOM and DINEOF in the southern area of the BoB HFR domain and that this regions are the same for the four experiments, large errors near the southern coastline.In addition of possible smoothing or noising, a discrepancies 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 throughout the BoB HFR domain.
Finally we quantify the difference between the LCS from the three methodologies and the REF-LCS by computing the spatial and temporal mean of the relative error defined in Eq. 3 and Eq.4 over all the LCS maps computed for the four experiments.
We summarize in Table 3 the results of these computations.As can be seen 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.Our results show that even when a large number of pixels is missing the FSLEs results still give an accurate picture of the oceanic transport properties and the velocity field filled by the three methods analyzed do not introduce artifacts in the Lagrangian computations. .

Residence times
Other Lagrangian quantity suitable to describe transport process is the residence times (RT) (Buffoni et al., 1996;Lipphardt et al., 2006;Hernández-Carrasco et al., 2013).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.
To compute RT particles are initially located on the grid points of the HFR velocity field and are integrated in time during 14 days.We only consider 14 days owing to the short period of the 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 pre-selected boundary after the 14 days of integration.
As done for the LCS, a first comparison is performed looking at the spatial distribution of RT obtained from the three different filled HFR.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, Figure 13 shows two maps of RT for the experiment A and C corresponding to April, 04, 2013, 2013 18:00. In    In order to reveal regions where the effect of the reconstruction on the RT is the strongest we compute the time average of the relative error for all the experiments (Fig. 15).This figure shows that, as seen for the LCS comparison, regions of accumulated large and small errors are the same for SOM and DINEOF and that this spatial distribution of the error is quite similar for the four experiments.OMA shows different regions with large errors, presenting the greatest number of points distributed through the HFR domain with large error.A common feature is that the three methods yield large errors near southern coastline and in the west and northwest area of the HFR domain.It suggests that the three methods are removing some dynamical feautures in this region.For points located in central regions (around 2 • W and 44 The four experiments based on different grouping of missing values demonstrate that even for spatially severe and persistent gaps, the LCS diagnostics obtained from the FSLE fields and the residence times applied to gap-filled high-resolution HFR surface currents are robust representations of surface dynamics in coastal basins.Our results show that even when a large number of pixels is missing the FSLEs computations still give an accurate picture of the oceanic transport properties and the velocity field filled by the three methods analyzed here are not introducing artifacts in the Lagrangian computations.The robustness of the Lagrangian diagnosis against the missing data and posterior filling process could explain the low relative error in the FSLE and LCS computations (Hernández-Carrasco et al., 2011).
While DINEOF presents the lowest errors in the Eulerian comparison of the velocity field, SOM is the method with lowest errors in the trajectories, LCS and residence times computations.DINEOF and SOM reconstruction methods are based on the patterns of the velocity extracted from statistical concepts.These method strongly depends on the number of modes/neurons used in the reconstruction.The larger the number of modes/neurons, the less smoothed the pattern will be.However even when using a large size of the neural network or number of modes, these methods are prone to filter the velocity field removing some small scale dynamical processes, and in some cases the resulting patterns are a smoothed representation of the real dynamics.DINEOF and SOM also depend on the choice of the modes or patterns.
In the SOM methodology the fact that each neuron is composed of velocity fields at three different times could increase the probability of using a more suitable pattern for the reconstruction, in particular, in the cases where there is a large number of missing points, i.e., experiment C (failure of antenna).This could explains the many points out of the main cloud in the scatter plot shown in Figure 5 for DINEOF for the experiment C.These outliers can propagate the error in the trajectories computations, explaining the difference between the Eulerian and Lagrangian comparison of the DINEOF and SOM methods.
On the other hand OMA is a geometrical approach that also depends on the choice of the modes and the number of modes (combination of irrotational and incompressible field configuration and the spatial shape of the HFR domain) used to project the velocity field.In general since only a finite number of modes are computed, we are introduced an arbitrariness in the selection of the modes, that is, the real velocity is projected in a subspace of modes that are either tangent to the coastline or in the open boundary.Using a small number of modes could affect the results of the reconstruction obtaining a velocity field far away from the HFR data but with very simple features.In the situations where the flow patterns are simple and quasi-permanent the use of a few modes is able to capture these features or major such dynamics.In contrast, if a large number of modes is used, the reconstructed field matches the HFRr data, but the data are not sufficiently filtered and smoothed, introducing some dynamical features or even artifacts, which could not be presented in the real flow.In our case, the coastline in the area of the HFR coverage is relatively large originating an increase in the number of modes tangent to the coast to the detriment of the OMA open boundary modes, therefore, conditioning the resulting inferred velocity field.As explained in section 3.1.the constraints applied in the OMA can limit the representation of localized small-scale features as well as flow structures near open boundaries.Besides, as discussed in Kaplan and Lekien (2007), difficulties may arise when dealing with gappy data, especially when the spatial gap size is larger than the minimal resolved length scale (that explains why the results are worst when there is just data from one antenna).The study of the individual sources of errors in the gap-filling methods is not trivial and that it depends on the coherence and persistence of the dynamical features as well as on the full knowledge of the dynamical scales involved in the tracers motion.Further analyses are needed to answer these questions, however this is out of scope of this manuscript and should be addressed in future studies.
In The time period selected here represents the dynamical conditions given in spring time (April).This is a transition period between the winter months (where the persistent Iberian Poleward Current dominates) and summer months where low energy and stable conditions are observed (L. et al., 2014).Applying the methodologies during the transition periods, characterized by high dynamical variability, allowed to avoid specific circulation patterns and ensure the results of this study are more easily generalized to other areas.However further analysis are required to corroborate this affirmation.
Concerning the temporal distribution of gaps, in this study, we have introduced gaps in the 50% of the time period which is much higher than the 10-20% of the accepted time failures in a well-functioning system.Nevertheless, this work is more focused in the spatial gaps and further analysis should be done in the future to analyze the limit of applicability of each method regarding the temporal horizon.
All the results reported in this study can not be extended to HFR systems working at different frequencies and resolutions.HFR working at different frequencies captures dynamical features at different scales that could be altered during the reconstruction process.Further studies using data from different HFR systems and regions characterized with different dynamical conditions should be performed to address these questions.For instance as said above some SOM and DINEOF modes used in the reconstruction are smoothed representations of the real dynamics and they could remove some small scale dynamical processes.However, these results are illustrative of the great performance of the gap-filling methodologies in providing reliable HFR velocities for a Lagrangian assessment of ocean coastal dynamics.
Contrary to SOM, and to a lesser degree to DINEOF, one of the main advantages of using OMA is that this method does not need a long time series of data and therefore it can be used immediately after installing the HFR system.Also, OMA allows the reconstruction of total velocities in areas of large GDOP, which also represents an advantage since it enables a larger spatial coverage.Besides DINEOF has the advantage of not requiring a long training dataset and do not having a subjective parameters.On the other hand, SOM is able to introduce nonlinear correlations in the computation of the patterns as well as to manage data sets with missing values.
These experiments also show that the developed SOM methodology is suitable to accurately reconstruct HFR velocity fields.
The results indicate that the SOM is able to encode properly the dynamical patterns present in turbulent flows.This is a very promising result and opens up new possibilities for applying this methodology for the inference of HFR at the dates with not available data.Although the performance of SOM and DINEOF methods is high, it is worthy to note that in this study we have used a simple algorithm based on the HFR dataset their own, which will be improved in future works performing the analysis on HFR velocities coupled with other oceanic variables.Moreover, a more rigorous sensitive analysis is needed to know the temporal coverage needed to obtain reliable reconstructed velocities field and to optimize the algorithm in terms of errors and computational cost.
A good approximation to obtain an optimal reconstruction of the velocity field could be the combination of this methodologies.For instance, firstly filling the gaps of the radial velocities by using SOM or DINEOF and then reconstructing the total velocities through OMA.It would allow to obtain a wider coverage without removing small scale features.
Two long-range CODAR Ocean Sensor SeaSonde HFR sites are operational in the SE BoB, since 2009.Both antennas emit at a 40-kHz bandwidth centered at 4.5 MHz frequency and they are part of the Basque in situ operational oceanography observational network owned by the Directorate of Emergency Attention and Meteorology of the Basque government's Security

Figure 1 .
Figure 1.Map showing the location of the two antennas (Matxitxako and Higer) of the BoB HFR systems.The grid points of the total HFR velocity field are plotted in blue.

For
all the cases data gaps are introduced in 50% of the timesteps.The proposed gap scenarios are: (A) Bearing gaps, generated by randomly distributing 10 of the elements of KMA group 6 (Figure 2); (B) Range decreasing generated by randomly distributing 10 of the elements of KMA group 4; (C) Antenna failure: in this group failures of one of the two antennas were simulated by randomly removing the data from one of the sites; (D) Random gaps: gaps were randomly introduced in the time series and randomly distributed in 30% of the spatial HFR domain.K-means analysis shows that the situation

Figure 2 .
Figure 2. Example of the four prototypes of KMA groups of gap distribution scenarios obtained by the KMA analysis applied to HFR availability-absence matrix for 2014.(1) representing a good data coverage scenario.(2) accounting of the failure of one of the two antennas.(3) characterizing a range decrease in both antennas.(4) representing bearing coverage decreasing in one of the antennas.Matxitxako radials are plotted in red and Higer radials are plotted in blue.(See figure S1 in supplementary material for the complete KMA results lattice).

Figure 3 .
Figure 3. Percent of spatial coverage for each time step and experiment.Higer data is shown at the upper panels and Matxitxako data at the lower panels in all cases.Black circles, percent of good data for the original fields.Red dots, percent of good data for each of the experiments.X-axis is the time step of the time series for April 2012, April 2013, April 2014.

Figure 4 .
Figure 4. Maps of temporal coverage (percent of the total number of hours in the 3-month period analysed) for Higer (upper panels) and Matxitxako (lower panels), for the REF series and all the experiments.
four experiments are shown inFigure (5)  andFigure (6)  respectively.In the case of OMA the cloud of points are, in general, widder with more points separated from the baseline.In general the behavior of SOM and DINEOF models are very similar in terms of RMS, BIAS and R 2 as depicted in Tables1 and 2for the u and for the v component of the velocity, respectively, although SOM tends to underestimate currents and DINEOF to overestimate.OMA is showing the worse reconstruction, specially for experiment C. It is worth to note 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 0.43, 0.38 and 0.57 for SOM, DINEOF and OMA, being the antenna failure the most difficult gap scenario to properly reconstruct.In particular, the scatter plot of original vs. reconstructed zonal velocities from DINEOF in the case of the Experiment C (Figure5second column, third row) shows many points out of the main cloud.It suggests that the gap scenario C

Figure 5 .
Figure 5. Scatter plot of the SOM (first column), DINEOF (second column) and OMA (third column) models vs observations for the four experiments (rows) for the zonal velocities.

Figure 6 .
Figure 6.Scatter plot of the SOM model (first column), DINEOF (second column) and OMA (third column) models vs observations for the 4 experiments (rows) for the meridional velocities.

Figure 7 .
Figure 7. Average of the separation distance between trajectories computed from the filled HFR using the three methodologies and reference HFR velocities as a function of time.Different limits in the vertical axis of the figure have been used for a better data representation.The error bar is the confidence interval (one standard deviation) of the spatial average of D(t).

Figure 8 .
Figure 8.Time 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.

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

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

Figure 12 .
Figure 12.Maps of the time average of the pixel by pixel absolute relative error of the LCS computed from the three filled with respect to LCS obtained from the reference velocity field.White color indicates points where LCS have nor been identified.
Time evolution of the spatial average over the HFR domain of the Absolute Relative Error, ARE, (Sec 4.1), of the attracting LCS obtained for each method and experiment with respect to the LCS computed from the reference velocity field is plotted in Fig. 11.A common characteristics for the three methods is the 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 mixing 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 in the middle of April 2013 and 2014.All the methods are more affected by the missing points given in the experiment C, in accordance with the Eulerian comparison.

Figure 13 .
Figure 13.Snapshots of RT computed from the three filled and the reference HFR for the experiment A and C at April 04, 2013 18:00..

Figure 14 .
Figure 14.Time series of the spatial average of the pixel by pixel absolute relative error of the three methodologies derived RT with respect to the REF-RT.

Figure 15 .
Figure 15.Maps of the time average of the pixel by pixel absolute relative error of the RT computed from the three filled with respect to RT obtained from the reference velocity field.
general, regions less affected by the missing points are located in the middle of the HFR domain.The large coverage of data can explain this spatial distribution.Moreover, large errors are concentrated in the west part of the HFR domain, likely induced because of the dynamics in this region is particularly more turbulent than in other areas over the SE BoB, increasing the effect of the velocity error on the trajectories of particles.The low number of gaps in the central part of the HFR domain in the experiment B and due to separation of the missing points randomly distributed in the experiment D can explain the similarity in the Lagrangian diagnosis obtained for these experiments and for the three gap-filling methodologies.With regard to the degree of the effect of the gap-filling methodologies on the reconstruction of the HFR signals depending on the type of gaps provided by the k-means gap patterns, experiment C, accounting for the failure of one of the two antenna, yields the largest errors in all the Lagrangian diagnosis according to the large Eulerian error.
Author contributions.I.H-C.conceived the idea of the study with the support of A.O, A.R and L.S; I.H-C developed the SOM methodology with the support of A.O; I.H-C produced the SOM HFR velocities and conduced all the Lagrangian calculations; A.O perform the Eulerian in supplementary material).Groups 1, 5, 8, 9, 10, 14 and 15 are considered as representatives of good data coverage; groups 4, 6, 7, 11 and 13 contain most of the common scenarios of spatial data-gap; and finally, groups 2, 3, 12 and 16 show situations with no data (or very few) for any of the two antennas (no totals currents can be produced).Examples of good coverage and gap scenarios corresponding to individual antenna failure, range reduction and bearing are shown in Figure2.
Since the goal of this work is to evaluate different gap-filling methodologies in real situations, the different groups representing observed gap types are used to introduce artificial gaps in April 2012, April 2013 and April 2014.Besides, randomly generated spatial gaps are also included to cover different failures that are not related to the typical situations but could occur.

Table 1 .
Root Mean Square Error (RMS), Mean Bias (MB) and Correlation Coefficient (R 2 ) of the zonal velocities (U) obtained using reconstructed HFR from the three models and for the four experiments respect zonal velocities derived from the Reference HFR currents

Table 2 .
Root Mean Square Error (RMS), Mean Bias (MB) and Correlation Coefficient (CC) of the meridional velocities (V) obtained using reconstructed HFR from the three models and for the four experiments respect zonal velocities derived from the Reference HFR currents Molcard et al. (2009)16)om the beginning of the trajectory simulation.DINEOF and SOM methods show similar results, with separation distances being around 1km after 24 hours of simulation and 3 km after 72 hours.Values of the distance for SOM are slightly better, in all of the experiments.On the contrary, this distance reachs values of 5 km after 24 hours in the case of OMA.The values of D obtained in our computations after 24 hours of simulation are smaller than those reported inBellomo et al. (2015);Kalampokis et al. (2016);Molcard et al. (2009)for separation distances between real and virtual 5 drifters trajectories.It is important to keep in mind that in these experiments the trajectories are computed from different source data (HFR and drifters) and sampling (spatial resolution) as well as different regions (Mediterranean

Table 3 .
Mean absolute relative error (ARE) and Mean bias (MB) of the LCS obtained using reconstructed HFR from the three methodologies with respect to the REF-LCS for the four experiments this figure we color the initial position of particles in the BoB attending to the time they need to scape the HFR domain.Regions from particles with short residence times are shown in cyan/blue, and initial positions of

Table 4 .
Mean absolute relative error (ARE) and Mean bias (MB) of the RT obtained using reconstructed HFR from the three methodologies with respect to the REF-RT for all the experiments It can be seen that except at the beginning of April 2012, in general, values of ARE for OMA-RT are larger than SOM and DINEOF-RT.In particular at the end of April 2012 and at the middle of April 2013 and April 2014.Again, experiment C turns out as the more critical situation of gap distribution to be reconstructed, reaching the largest values of ARE close to 0.25.
(Yaremchuk and Sentchev, 2009;Fredj et al., 2016) HFR velocities approximate better to 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 Table4the results of the computations of the ARE and ME by making averages over all the RT values derived from the three filled HFR velocities to be validated with the REF-RT.A total of 156240 RT values are used in the computations.Although the absolute relative error, is slightly larger in the case of the OMA-REF than using SOM and DINEOF, these errors are low in all cases, being 0.16 the maximum ARE obtained among all the experiments.The absolute error, is smaller for SOM and DINEOF cases (between 0.12 and 0.13 for both cases) and slightly larger for OMA (between 0.13-0.16).Furthermore the ARE for SOM-RT and DINEOF-RT are quite identical.In all cases the MB are positive, meaning that RT are overestimated.Largest values of MB are found, in average, for OMA-RT (MB ranging from 0.7 hours for experiment D to 8.1 hours for experiment A, following by SOM-RT (MB ranging from 5.0 hours for experiment B to 7.8 hours for experiment C) and then by DINEOF-RT (MB ranging from 4.3 hours for experiment B to 5.4 hours for experiment A).We have investigated the reliability of the HFR gap-filling methodologies by studying the sensitivity of some basic Lagrangian metrics for the assessment of the dynamical properties when they are computed using reconstructed HFR velocities.The sensitivity test is carried out through four experiments including the most common scenarios of data gaps observed in HFR systems due to different environmental or system's own causes, i.e., hardware failures, communication interruptions, particular environmental conditions affecting the detection of the signal, sea state, etc.Due to the importance of Lagrangian properties of the ocean for coastal management, and in contrast to other comparative studies of different gap-filling methodologies that are based on Eulerian differences of the filled and original velocities(Yaremchuk and Sentchev, 2009;Fredj et al., 2016), we have focused this comparison on the Lagrangian diagnosis.We have seen that the relative errors found for the reconstructed velocity field with respect to the original HFR field, (up to 43%, 38% and 57% in the case of SOM, DINEOF and OMA, respectively, for the worst gap-scenario corresponding to one antenna failure) are reduced in the Lagrangian diagnosis.The largest errors found between the four experiments and for the LCS computations are up to 22%, 22% and 27% and up to 13%, 13% and 16% for RT, when using SOM, DINEOF and OMA reconstructed velocities, respectively.