Effect of small scale transport processes on phytoplankton distribution in coastal seas

Coastal ocean ecosystems are major contributors to the global biogeochemical cycles and biological productivity. Physical factors induced by the turbulent flow play a crucial role in regulating marine ecosystems. However, while large-scale open-ocean dynamics is well described by geostrophy, the role of multiscale transport processes in coastal regions is still poorly understood due to the lack of continuous high-resolution observations. Here, the influence of small-scale dynamics (O(3.5–25) km, i.e. spanning upper submesoscale and mesoscale processes) on surface phytoplankton derived from satellite chlorophyll-a (Chl-a) is studied using Lagrangian metrics computed from High-Frequency Radar currents. The combination of complementary Lagrangian diagnostics, including the Lagrangian divergence along fluid trajectories, provides an improved description of the 3D flow geometry which facilitates the interpretation of two non-exclusive physical mechanisms affecting phytoplankton dynamics and patchiness. Attracting small-scale fronts, unveiled by backwards Lagrangian Coherent Structures, are associated to negative divergence where particles and Chl-a standing stocks cluster. Filaments of positive divergence, representing large accumulated upward vertical velocities and suggesting accrued injection of subsurface nutrients, match areas with large Chl-a concentrations. Our findings demonstrate that an accurate characterization of small-scale transport processes is necessary to comprehend bio-physical interactions in coastal seas.


Results
Relevance of small scale processes to properly map coastal transport. We first provide evidences to support the reliability of the HFR currents for our computations. The Lagrangian validation has been performed using data from 8 drifters trajectories available in the domain of interest (HFR area of coverage) during the period of study (July 2012-July 2014). We compare synthetic drifter trajectories integrated in the HFR velocity field and real drifter trajectories by computing the distance (D) between each virtual drifter trajectory and the real one 51 . The largest values of D were found for drifters deployed in October 2012 and April 2013 reaching a separation distance of 7-8 km after 24 hours (see Fig. S2 a in Supplementary Information). The smallest D values (~1 km) were found for the drifters deployed in September 2012. This difference in the D values could show the different dynamical conditions between both months. Nevertheless, the mean distance of separation, i.e. averaging over all the trajectories, after 24 hours is about 4 km (see Fig. S2b). It demonstrates a good accuracy of the HF Radar measurements as well as of the particle trajectories derived from it.
The performance of HFR-LCS (Lyapunov ridges extracted from HFR velocities, -see Materials and Methods) in structuring the flow is shown by advecting in the HFR velocity field two sets of virtual neutrally buoyant particles initially deployed on the northern and southern flank of a given LCS that lays across the area from east to west in October 28, 00:00 UTC (Fig. 1a). Although the location and magnitude of this LCS evolve in time, the LCS persists for several hours manifesting the presence of a coherent transport barrier preventing both sets of particles to be mixed up. A meridional LCS is formed and maintained during the simulated period, limiting water exchanges between the coast and the open ocean ( Fig. 1a-h). The significance of these HFR-LCS in organizing the transport in this coastal region is also supported by comparing their evolutions with real drifters trajectories. It constitutes an independent and observed dataset to further validate our approach. In Fig. 1, the position of a real drifter is superimposed on the backward FSLE field at every snapshot. The drifter is initially located close to an attractive LCS and it remains in the vicinity of the structure during the entire period. The same validation has been performed for other 7 drifters with similar results (see Movies in Supporting Information). To evaluate how backward LCSs constrain the drifters motions we have compared the probability distribution function (PDF) of the values of the backward HFR-FSLEs at the position of the real drifter with the PDF of the FSLE values for the entire HFR domain (Fig. S1c). As compared with the peak at low FSLEs found in the rest of the domain, a clear peak of the PDF at high FSLE is found at the drifters locations. This larger occurrence of high FSLE at the drifter positions shows that drifters are attracted to backward LCSs. Note that in contrast to other authors 51,52 we use Lagrangian averaged metrics, i.e. LCSs, instead of trajectories to validate the HFR velocity field and the derived Lagrangian computations.
Both tests based on virtual particles and using real drifters demonstrate that even for restricted domains (~thousands of km 2 ) and small temporal scales (few hours to a few days), the LCS diagnostics obtained from the FSLE computation applied to high-resolution HFR surface currents are robust representations of surface dynamics in coastal basins.
Following other authors 17,45,53 we now evaluate the dynamical importance of the scales captured by the HFR by computing particle dispersion at different spatial scales using Eq. 1 (see Material and Methods section). The averaged FSLE curve (Fig. 2) shows that λ(δ) is not constant over the range of scales covered by the HFR (scale dependent). The best-fitting of the non-averaged (or individual) FSLE curves returned slopes spanning −0.64 to −0.73; the best fit of the averaged FSLE curve shows a slope of −0.68 with a correlation coefficient, R 2 , of 0.97. This indicates that the average scaling law of the relative dispersion over the scales ranging from 3.5 km to 25 km is associated to a Richardson turbulent diffusion 54,55 . It means that the main contributors to the separation rate are structures with size comparable with the separation itself, therefore the dispersion regime is local.
This shows that the lateral relative dispersion, at the surface, is governed locally by small-scale structures and not only by meso-and larger-scale structures. It also suggests that HFR is partially resolving submesoscale i.e. it is capturing the largest (~3.5-10 km) submesoscale features of the flow. Targeted fitting analysis on specific scales shows that over 3.5-8 km the slope is smaller (−0.60) than for scales ranging 8-25 km (−0.71). The dependence of λ(δ) at δ larger than 25 km cannot be explored due to the limited area covered by the HFR. Note that those scales (<25 km) are not resolved by the geostrophic velocities inferred from altimetric Sea Surface Height (SSH), thus the use of HFR velocity fields is necessary to accurately study transport processes in this region. From the fitting in the Richardson regime (λ(δ) ~ε δ ⋅ − 1/3 2 /3 ) 55 we obtain a mean turbulent dissipation rate ε of 2.5 ⋅ − 10 m /s 9 2 3 , which is similar to values reported in most oceanic regions by Corrado et al. 17 , using drifters.
The slope of the FSLE curve for synthetic drifter trajectories integrated with HFR velocity field is consistent with previous observational and modeling studies, with a slope similar to those derived from surface drifters 16,26,56 and from numerical simulations 22 . It further reinforces the fact that the velocity field of the HFR is adequate to study dynamics over the upper submesoscale (3-10 km) and the lower mesoscale (10-30 km).
The λ obtained in our computations is one order of magnitude smaller than the one of Haza et al. 45 , (λ = 4 to 7 − days 1 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 [3.5-25 km] 17 . Comparing with other observational studies in the Mediterranean Sea and at this range of scales we found values of the FSLE curves of the same order (O(1) days −1 ). Schroeder et al. 56 , using drifters deployed in the Ligurian Sea reported slightly larger values of λ (between 0.9 and 0.4 days −1 ) at spatial scales between 3 and 30 km following a Richardson slope. In the Adriatic Sea similar values of λ (0.7-0.2 days −1 ) were reported for the scales range studied here 26 . Other studies based on drifters in the Gulf Stream 15 have found higher values of λ(δ). In addition of possible observational biases, various regional discrepancies could explain these different values of FSLE. To name just a few, a strong coastal jet, alternating tidal currents, prominent topographic boundaries, specific wind regimes, etc… will influence transport dynamics and/or affect near-surface measurements. Impact of lateral stirring on coastal Chl-a distribution. Next we characterize the biological response to transport processes by comparing the HFR-LCSs with patterns of satellite Chlorophyll-a (Chl-a) at 1 km resolution as a proxy of surface phytoplankton concentrations. Figure 3 shows four snapshots of the spatial distribution of surface Chl-a together with the HFR-LCSs (in black) corresponding to winter and early spring of 2013 ( Fig. 3a and b) and 2014 ( Fig. 3c and d). Spatial distribution of Chl-a is very heterogeneous within this relatively small area, exhibiting important phytoplankton patchiness. As observed in each snapshot the attractive LCSs shape the spatial distribution of Chl-a. High concentrations of Chl-a are constrained ( Fig. 3a and b) and stirred ( Fig. 3c and d) by Lyapunov lines.
Phytoplankton blooms are commonly observed in this region during winter and early spring 33,57 . Figure 4a shows the temporal evolution of the daily Chl-a concentration spatially averaged over the whole HFR area (magenta annotation, Fig. 1a). We also display the time-series averaged over small northwestern and southeastern boxes (magenta boxes, Fig. 3a) to study the spatial propagation of the bloom. Boxes are chosen following the main circulation pattern (northwestward) given by the first mode of the EOF analysis of the HFR currents 29 .  As observed in the southern control box, the maximum of Chl-a appears at the end of January 2013 and for the northern control box at the end of February -early March. A less intense bloom is also present during 2014 with a slight shift of the maximum between the two control areas. A possible scenario is that nutrients-enriched coastal waters are advected by the northward flow of Atlantic waters, intruding the HFR area and mixing (small-scale structures) with the nutrient depleted Liguro-Provencal-Catalan current.
As a measure of mixing activity, we compute spatial means of FSLE values for the whole region and for both control areas (Fig. 4b). FSLE maxima coincide with maxima of Chl-a. Large FSLE values appear first in the south box and then in the northern one, perfectly matching the meridional order of appearance of Chl-a maxima. The best non-linear fitting between daily Chl-a and FSLE values over 2 years gives a correlation coefficient of 0.70 (Fig. 4e) for a polynomial fitting of second order. High mixing activity produced by the stirring of surface currents is involved in the generation of water masses with large values of phytoplankton. This has been also observed by Calil and Richards 11 , where large distribution of phytoplankton in oligotrophic areas was related to the dynamical mesoscale structures inferred from altimetry. By contrast, other studies showed that large surface mixing activity has a negative effect on phytoplankton productivity in upwelling regions 12,42 , which was explained as a "regional-scale" lateral-induced loss of nutrients 14 . This strong positive correlation between FSLE and Chl-a is an indicator of the tight bio-physical interactions at the small-scales revealed by HFR. Other factors could be involved in the biomass enhancements, namely, the seasonal development of the mixed layer depth (MLD) and of the thermocline. Lavigne et al. 34 , based on observations have reported that this region exhibits a very shallow MLD in summer, reducing the supply of nutrients from deep waters, and a deeper MLD in winter, as a consequence of the specific atmospheric conditions, favoring the fertilization of surface layers. Further comparison of the MLD with the nutricline and euphotic depth has shown that primary production was limited by nutrient fueling from below rather than by light availability 34 . This suggests that the growth of phytoplankton is controlled by either a MLD deeper than the nutricline, or by the vertical advective transport of nutrients induced by the small scale processes within the mixed layer. Indeed, previous works have shown that the supply of nutrients to the euphotic layer is related to the submeososcale processes, i.e., the stretching of filaments through horizontal advection (see Mahadevan 7 for a review). While our results clearly suggest that, in this coastal region, high mixing activity is associated with abundant phytoplankton, the precise mechanisms remain unknown. We also computed the FSLE using geostrophic currents from Altimetry (SSH-FSLE) (see Materials and Methods) and we found that the temporal evolution of spatially averaged SSH-FSLE did not capture the seasonal change of mixing activity seen by the HFR (Fig. 4c). No significant correlation is found between the averaged values of geostrophic FSLE and Chl-a over the whole region nor the averaged values of both control areas. Accumulated effect of small scale sea surface distortion on phytoplankton patchiness. To further study the mechanisms explaining the relationship between phytoplankton patchiness and transport processes, we now exploit a geometrical quantity related to the contraction and expansion of the flow. If a three-dimensional flow is incompressible then, following the continuity equation, the three-dimensional divergence of the velocity field is zero and one can use the two-dimensional horizontal divergence at the surface to differentiate areas of contraction and expansion of the fluid flow. If the ocean flow is strictly two-dimensional at the sea surface, divergence will be zero and the corresponding ridges of the FSLE field correspond to lines of stretching and folding. However, in truly 3D flows, when vertical motions exist, the horizontal divergence will be non-zero and the horizontal convergence/divergence zones at the surface are related to stretching and folding areas as well as to vertical motions. Positive divergence indicates upward motion and negative divergence is associated with downward motion. These conclusions do not hold for altimetry-derived velocities because the geostrophic approximation results in horizontal currents (neglected vertical velocities) of zero-divergence.
Assuming that the divergence of surface velocity field derived from HFR data is non-zero, we analyze the accumulated effect of convergence and divergence of a fluid parcel by computing the Finite-Domain Lagrangian Divergence (FDLD) given by Eq. 5 (see Materials and Methods). The integration duration is chosen according to the typical time-scales required for phytoplankton communities to respond to a new supply of nutrients, which is around 5 days in frontal systems 58 . Snapshots of FDLD for two different days displayed in Fig. 5(a and b) show complex patterns, attesting of the high spatial variability of convergence and divergence zones.
Comparing spatial patterns in Figs 3d and 5b, we find that convergence lines given by minimum (negative) values of FDLD (ravines in FDLD field) identify attractive LCSs, thus revealing dynamical structures that affect transport. Figure S3 in Supporting Information shows further exemplary comparisons between both Lagrangian structures. Figure 4i reveals a linear fit as the best fit found between both dynamical (backward FSLE) and geometrical (FDLD) spatially-averaged quantities with a negative correlation (R 2 = −0.79). Negative values of FDLD suppress expansion and enhance trapping, giving rise to regions of clustering. In those regions, high values of Chl-a are concentrated along the negative lines of FDLD (Figs 3d and 5b) which match the attractive LCS. In this case also a polynomial fit of second order was found as the best fit between Chl-a and negative FDLD with a correlation of R 2 = 0.66 (Fig. 4g). Note that since FSLE ridges also identify regions of high shear, only the well-matched patterns of FDLD and FSLE can be used to distinguish local contraction and expansion motions. The numerical study by Jacobs et al. 50 , from models at different resolutions reveals that surface material clustering is initially dominated by these small scale features. We have seen here how the clustering effect of small scale processes contained in HFR currents can affect horizontal phytoplankton distribution.
Next we explore the effect of vertical motions associated to the dynamical structures derived from LCS and FDLD on phytoplankton distribution. Calil and Richards 11 showed that regions of pronounced gradients of relative vorticity, related to backward FSLEs, are expected to exhibit large vertical velocities. The scatterplot displayed in Fig. 4h shows good correlations between spatial averaged forward FSLE and positive FDLD (R 2 = 0.78). The effect of these geometrically-induced vertical motions on phytoplankton concentrations is shown in Fig. 5a. We find that regions of large Chl-a concentrations (Fig. 3c) are associated to regions of positive FDLD values (Fig. 5a) that indicate accumulated upward vertical velocities. Two distinct mechanisms could explain the fact that extreme FDLD values (both positive and negative) are associated with phytoplankton patches. Positive FDLD extrema indicate divergence and the predominance of vertical upward motions of subsurface nutrients, thus promoting locally phytoplankton growth. Negative FDLD extrema mean convergence and horizontal clustering, thus accumulating surrounding phytoplankton. To disentangle vertical and horizontal dynamics and to determine the preponderant mechanism, our results advocate for carefully analyzing the sign of the FDLD extrema as well as the relationship between FSLE and FDLD fields.
Time-series of the spatial averages of FDLD (Fig. 4d) show that maxima of FDLD associated to high positive divergence (upward motions) occur concomitantly with the growth of phytoplankton (Fig. 4a) and with the peak of mixing activity (Fig. 4b). This notable seasonal variability of FDLD is not seen for Eulerian (instantaneous) divergence (Fig. 4d). Positive Lagrangian divergence prevails over negative values during the phytoplankton bloom periods. The scatterplot of Chl-a vs. FDLD displayed in Fig. 4f shows the best non-linear fitting (polynomial fit of second order) with a positive correlation between both quantities (R 2 = 0.68).
Note however that this significant correlation between Chl-a and FDLD does not hold when the comparison is made pixel-by-pixel. This misfit could be due to the fact that those observations (currents and Chl-a) are not exactly concomitant as they originate from different remote sensors or to the tentative comparison between an Eulerian (Chl-a) and a Lagrangian (FDLD) quantity. Also daily Chl-a concentrations derived from satellite satellite suffer from gaps due to cloudiness (thus reducing the number of pixels available for the comparison). To circumvent this issue, we adopted a novel methodology inspired from the Lagrangian framework: we integrated the Chl-a scalar field along the trajectories provided by the HFR during the same integration period of the FDLD to estimate the "cumulative amount of positive FDLD" seen by each patch of Chl-a along its drift. The correspondence between both variables is improved when comparing two snapshots of accumulated Chl-a ( Fig. 5c  and d) with the corresponding FDLD patterns ( Fig. 5a and b, respectively). Regions where the accumulated Chl-a is large match well regions of high negative FDLD values (where clustering is strong) as well as high positive FDLD values (where upward motions are predominant). To quantify the relationship between vertical motions and phytoplankton concentrations we compute the correlation coefficient between values of accumulated Chl-a and values of positive FDLD. We obtain a value of R 2 = 0.56, which is much larger than when using the Eulerian Chl-a (R 2 = 0.25).

General Discussion and Conclusion
We investigate coastal transport using two complementary Lagrangian diagnostics, FSLE and FDLD, both computed on a high-resolution velocity field derived from HFR. The combination of both quantities allows us to distinguish regions of high shear, convergent regions where clustering dominates, and areas with persistent vertical motions. These analyses are then used to characterize the physical processes which impact biological activity and patchiness in coastal regions.
The scale-dependence of HFR-FSLE provides an evaluation of the size of the dynamical features captured by the HFR affecting the transport properties. The slope of the averaged FSLE curve over the scales comprised between 3.5 km and 25 km is associated to a Richardson-like dispersion. This slope is in agreement with previous observational studies with drifters 16,26,56 . At the range of scales measured by HFR, surface relative dispersion is locally controlled by structures in the upper submesoscale range (~3.5-10 km) and in the lower mesoscale (10km-25km).
HFR-LCSs are prominent structures that organize transport and dispersion of coastal waters in the Ibiza channel, impacting physical and biological coupled processes at regional scales. HFR-LCSs constrain the motion of virtual particles and real drifters, and shape the distribution of Chl-a concentration. Lehahn et al. 9 , showed that mesoscale fronts obtained from SSH-FSLE separate regions of different Chl-a concentrations in the open ocean. However, in coastal basins, SSH-FSLE is not able to capture the seasonality of the mixing activity nor the spatial distribution of LCSs provided by the HFR-FSLEs. This reveals that small scale processes are important to understand phytoplankton patchiness in coastal basins. The influence of small-scale coastal dynamics on surface phytoplankton distribution in coastal areas can be properly studied considering (i) small spatial-scales (  0 0345 . vs. . 0 125  resolutions of HFR and SSH, respectively), (ii) high frequency data (hourly vs. daily velocities of HFR and SSH, respectively) and (iii) the total velocities (geostrophic + wind-driven + Stokes drift) currents of HFR data vs. geostrophic component of the SSH velocities). As an additional example, we compare the SSH-LCSs with the HFR-LCSs (see Fig. S4a and S4b in SI) at the same spatial resolution of 0.4 km. Although both computations identify an eddy-like structure, the one derived from altimetry is shifted westward as compared to the eddy detected by HFR. Similarly, the finer filamentary structure provided by HFR-FSLE is not well described by SSH-FSLE.
Using FDLD, we show that such small scale structures can be further distinguished in those that likely produce clustering and those that are rather associated with large accumulated upward velocities. On one hand, regions of highly negative values of FDLD are convergent and prone to clustering; thus they accumulate the surrounding phytoplankton standing stocks. On the other hand, regions of positive FDLD are associated to the prevalence of upward vertical velocities that promote new production of phytoplankton due to enhanced nutrient supply. By comparing the fields of our dynamical Lagrangian diagnostics with high resolution maps of Chl-a, we found that regions of high Chl-a match regions with both positive and negative extrema of Lagrangian divergence, while only the latter coincide also with attractive LCSs. In our area of study and over typical time-scales of a few days, positive divergence predominates during the periods of phytoplankton blooms in fronts, suggesting that accumulated upward velocities produced by the highly divergent horizontal flow favor phytoplankton growth through the upwelling of subsurface nutrients, thus controlling the biological dynamics.
The relationship between the physical (FSLE and FDLD) and the biological (Chl-a) variables is first evaluated through the significant correlation coefficients found by comparing spatially averaged values. The best fitting analysis between both variables shows good correlations for polynomial fitting of second order. This suggests that stirring (FSLE) and flow divergence (FDLD) have significant effects on phytoplankton distribution when they reach extreme values; i.e. when mixing activity and divergence/convergence are high, only high values of Chl-a are observed. When the mixing activity or the divergence/convergence are low to moderate, their effects on phytoplankton is not so clear. The relationship between these two physical drivers and Chl-a is therefore symmetric; they can be viewed as enhancing factors for planktonic biomasses.
Furthermore we explored the spatial agreement of the FDLD with Chl-a estimating the Lagrangian accumulation of Chl-a along the trajectories of the fluid parcels. We go to the accumulation of Chl-a because we want to see the accumulative process of nutrients injection. This approach returned a much better correlation between both variables than using the "Eulerian" Chl-a. Of course a tight absolute relationship between Chl-a and FDLD is not expected since many other key biological processes for phytoplankton (e.g. heterogeneous light field, recycled production, zooplankton grazing, viral infections, etc…) were not taken into account here.
As ocean surface currents forecasting systems are being developed 5,59 , our methodology and the diagnostics exploited here could be applied on current predictions and thus serve to project future spatial distribution and magnitude of Chl-a in the coastal ocean. More generally, our study proves that the combination of both LCSs and FDLD computed from HFR is a powerful framework to study the effect of transport on biological quantities in coastal seas, as well as to localize convergence/divergence zones which are relevant for the tracking of debris accumulation or jellyfish aggregations, and for other coastal management activities, such as search and rescue missions and oil spill management.

HF Radar measurements and drifter data. A coastal HFR network combining two CODAR antennas in
the Ibiza Channel (IC) area provides hourly velocities up to 70 km from the coast of Ibiza with a radial resolution of 1.6 km (data available at http://www.socib.es). The coverage area of the HFR is shown by the magenta circle plotted in Fig. 1a. Cartesian velocities are provided for the period of June, 1st 2012 to July, 31st 2014 in a regular grid of 393 points over an area of approximately 3500km 2 with a resolution of 3 × 3 km and with an accuracy of ±0.5 m/s 29 . In order to be able to compute particle trajectories the gaps of the HF Radar velocity field have been filled using the OMA (Open-boundary Modal Analysis) algorithm developed by Kaplan and Lekien 60 .
Drifter data used for the validation of the HF-Radar velocities and HFR-LCS are part of different experiments performed by ICTS SOCIB (http://www.socib.es) in the Ibiza Chanel. Drifters position were provided through GPS positioning transmitted via GSM every 5 minutes.
Satellite data. We use daily AVISO Sea Level Anomalies at 1/8° of spatial resolution derived from a Ssalto/ Duacs multimission altimeter regional product released in 2016 to compute the Absolute geostrophic surface velocities. This product can be downloaded from (http://www.aviso.altimetry.fr/en/data/products/ sea-surface-height-products.html). It consists of surface currents computed using time-variable Mediterranean Absolute Dynamic Topography (η) calculated from mapped altimetric sea level anomalies combined with a mean dynamic topography computed for the Mediterranean Sea 61 . Then the zonal and meridional components of the Absolute geostrophic currents, ( φ θ U , V ), are obtained from the geostrophic relationship. The velocity field used here covers a period from June 2012 to June 2014.
We use ocean surface Chlorophyll data obtained from the ESA-CCI Remote Sensing Reflectance spectrum using a specific regional processing algorithm for the Mediterranean (MedOC4). The reflectance is the result of merging MODIS-Aqua, SeaWiFS and MERIS sensors. The product is distributed by CMEMS and can be downloaded from http://marine.copernicus.eu. This product measures the average chlorophyll content over the optical depth. Gridded daily data were used with a spatial resolution of approximately 1 × 1 km.

Finite Size Lyapunov Exponents and Coastal Lagrangian Coherent Structures.
We quantify horizontal transport processes by the Lagrangian technique of the Finite Size Lyapunov Exponents (FSLEs) 36,62 . FSLE is specially suited to characterize flows with multiple spatio-temporal scales 53,63 , and to study the stretching and contraction properties of transport in geophysical data 38,64 .
FSLE was originally introduced in the dynamical system theory to characterize the growth of non-infinitesimal perturbations in turbulence 36 by the following equation, λ δ τ δ δ = 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.
The explicit dependence of this expression of average FSLE on the separation scale δ allows to isolate the contribution of different spatial scales to particle-pair separation 53,63 , providing a characterization of the relative dispersion processes in terms of scaling laws. In contrast to relative dispersion metrics treating the time as the independent variable, λ(δ) (spatial scale, rδ, as independent variable) is not affected by the interference of pair-particles at different dispersion regimes. The value of r has to be small (r ¡ 2), in order to capture the relative dispersion from small coherent features and avoid the contribution of different scales of motion in the average, and not close to 1 to avoid aliasing problems related to the time step of particle advection at small scales 63,65 .
We compute the averaged FSLE with a total of 648 particle pairs trajectories released in the area of HF Radar velocity field for each δ ranging from 3.5 km to 25 km, and r = 2 . We repeat the same computations for 480 different initializations (homogeneously distributed through the two years of data). Particles escaping the HF Radar domain are not taken into account in the average computation. Only original pairs have been taken into account for the average. The particle trajectories in a 2D HFR velocity field are governed by, are the zonal and meridional velocities measured at coordinates = x y r ( , ) by the HFR. FSLE can also be used to unveil dynamical structures that act as transport barriers [37][38][39] . The calculation of the LCSs from FSLEs goes through computing the minimum time, τ, required for two fluid particles initially centered in r and separated by a distance δ 0 to reach a fixed final separation distance δ f . At position r and time t, the FSLE, denoted by λ, is given by: t l n r ( , , , ) ( In this definition of FSLE 0 δ δ = and r f δ δ = using a large value of r to adequately distinguish regions of extrema of the FSLE field. Note that, to have an explicit space-time dependence, averages are not performed here, in contrast with the original definition 1. The largest Lyapunov values concentrate along characteristic lines, Lyapunov lines, which can approximate manifolds of relevant hyperbolic points, the so-called Lagrangian Coherent Structures (LCS) 35,67 . Fronts, eddies and filamentary barriers to transport can be identified with these manifolds. Since LCS cannot be crossed by particle trajectories, such lines strongly constrain and determine fluid motion, organizing ocean transport.
To estimate the LCS from FSLE computations we use the algorithm described in Hernandez-Carrasco et al. 39 , adapted to finite domains, in this case the area of HFR coverage. The minimum time τ is computed by integrating the trajectories of the four neighboring points of the analyzed one located at r and by selecting the associated particle that separates faster δ f . The grid points whose four neighboring trajectories escape the domain before separating δ f are not taken in the analysis. The relation between δ 0 and δ f has been calibrated through previous experiments in order to obtain clear and robust LCSs in the relatively small domain of the HFR coverage. The total time of integration is chosen to be greater than the maximum value of the average residence time obtained in this area. In our computations the residence time is always less than 15 days so we used a maximum of 15 days for integrating the particle trajectories. We compute the FSLEs for all points x in a lattice of initial conditions, defined by a grid spacing coincident with δ 0 .
Numerically, we integrate Eq. (2) using a standard, fourth-order Runge-Kutta scheme, with an integration time step dt 10 = minutes. Since information is only provided in a discrete space-time grid, spatio-temporal interpolation of the velocity data is performed by bi-linear interpolation. FSLEs are obtained for the points r of a lattice with spacing δ 0 = δ. Initial conditions for which the prescribed final separation δ f = rδ has not been reached even when using all available times in the dataset are assigned a value λ = 0. Assuming that HFR velocity field is divergent, we accumulate the instantaneous horizontal divergence along a trajectory s x y t ( , , ) 0 0 0 in the limited domain, which we refer to as Finite-Domain Lagrangian Divergence, integrating and averaging its value over time as followed: x t y t t dt v ( , , , ) 1 ( ( ), ( ), ) , = − , is the time interval of the trajectory integration. As for the FSLE, we define the Finite-Domain Lagrangian Divergence (FDLD) as the integrated Eulerian divergence of the fluid particles over their trajectories within the finite area covered by the HFR. Note that in our computations t f may be different for each particle. All particles are integrated during a fixed time interval T, whose value must be chosen depending on the process studied and on the available domain (covered by the velocity field). To prevent trajectories to escape area during the integration, T must be lower or equal to the residence time (RT) in the region of interest. In our case, we use T = 5 days (typical response time of phytoplankton bloom in coastal fronts) when possible and T RT = otherwise, with RT ranging from 3 to 8 days in our region of interest. This Lagrangian measure has recently been applied in atmospheric flows by Tang et al. 47 : based on LiDAR observations, they compared the LCSs against the Lagrangian Divergence ridges, as a proxy of divergent and convergent regions, to finally interpret the 3D structures of the flow. We also refer the readers to Huntley et al. 49 , who performed an extensive comparison between Finite-Time Lyapunov Exponents (FTLE) and Lagrangian Divergence using modeled data. They proved that FTLE can be expressed as the sum of the stretch and dilation rates, with the latter being the averaged divergence along the trajectories of the fluid particles. A rough approximation of the relationship between classical Lyapunov Exponents and Lagrangian divergence is given in Section S3 in the Supplementary Information following Falkovich et al. 68 . Here we use HFR velocities, which represent the total velocity at the sea surface, and we compare these geometrical quantities with biological observations. In the case of the currents obtained from satellite altimetry the geostrophic approximation has to be considered. The geostrophic balance between horizontal pressure gradients in the ocean and the Coriolis force results in horizontal currents (zero vertical velocities) with a neglected horizontal divergence. Therefore the divergent/convergent effects cannot be observed by satellite altimetry.