Spreading of Lagrangian Particles in the Black Sea: A Comparison between Drifters and a High-Resolution Ocean Model

A Comparison between Drifters and a High-Resolution Model. Abstract: The Lagrangian dispersion statistics of the Black Sea are estimated using satellite-tracked drifters, satellite altimeter data and a high-resolution ocean model. Comparison between the in-situ measurements and the model reveals good agreement in terms of the surface dispersion. The mean sub-basin coherent structures and currents of the Black Sea are well reproduced by the model. Seasonal variability of the dispersion in the upper (15 m), intermediate (150 m) and deep (750 m) layers are discussed with a special focus of the role of sub-basin scale structures and currents on the turbulent dispersion regimes. In terms of the surface relative dispersion, the results show the presence of the three known turbulent exponential, Richardson and diffusive-like regimes. The nonlocal exponential regime is only detected by the model for scales <10 km, while the local Richardson regime occurs between 10 and 100 km in all cases due to the presence of an inverse energy cascade range, and the diffusive-like regime is well detected for the largest distance by drifters (100–300 km) in winter/spring. Regarding the surface absolute dispersion, it reﬂects the occurrence of both quasi-ballistic and random-walk regimes at small and large times, respectively, while the two anomalous hyperbolic (5/4) and elliptic (5/3) regimes, which are related to the topology of the Black Sea, are detected at intermediate times. At depth, the signatures of the relative and absolute dispersion regimes shown in the surface layer are still valid in most cases. The absolute dispersion is anisotropic; the zonal component grows faster than the meridional component in any scenario.


Introduction
The Black Sea (BS) is a semi-enclosed dilution basin with limited water mass exchange with the Mediterranean Sea. In its southwestern region, the BS is connected with the Mediterranean Sea by the Bosphorus Strait, through which a small amount of salty water flows into the BS. Despite its small size, the Bosphorus Strait plays a critical role for the hydrodynamics and climatology of the BS [1]. In the northwestern BS area, the inflow of freshwaters (low-salinity waters) from large rivers (e.g., Danube, Dnepr and Dnestr rivers) leads to the generation of a pycnocline near 100-150 m, which inhibits exchanges between the surface and deep waters [2]. This freshwater is rich in nutrients and is also contaminated with industrial and mining wastes. Major contaminants include oil residues, pesticides, hydrocarbons, nutrients, heavy metals and accidental oil spills. These contaminants pose a very high level of risk for the marine environment and coastal areas in the BS [3,4], and hence, it is critical to understand the dispersion properties of particles in this area. The Sub-mesoscales and mesoscales features, filaments, currents and fronts are responsible for physical transport in the ocean, whose understanding is crucial for a number of important applications [17,18]. The turbulent dispersion is affected by the presence of these features that induce a discrepancy between the measured, modeled and the theoretical (2-D turbulence theory for isotropic and homogeneous flows) dispersion laws. To provide an extensive discussion of this discrepancy, it is necessary to compare the theory with in-situ data and simulations. The study of turbulent dispersion in the BS and its contribution to transport processes is very important for many reasons. For instance, the BS is an enclosed Sea with limited water exchange with the open sea. Therefore, it is one of the most isolated seas characterized by large continental shelf where there is important biological production (nutrients, fish eggs, jellyfish). Another reason for studying the dispersion is to address priorities such as threats, oil and chemical wastes in this area, which are known to be present in the most contaminated regions.
Previous studies have been focused on the dispersion analysis in the world ocean in order to demonstrate the dispersion properties and to help to describe/predict the spreading of pollutants. Examples of this focused on the Gulf of Mexico [19], the Gulf of Stream [20] and the Mediterranean Sea [21][22][23][24].
Turbulent dispersion of floating particles is influenced by sea waves [25,26], tides [27] and circulation [19][20][21][22]. In this study, the role of the upper, intermediate and deep BS coherent structures on the water mass dispersion is investigated.
The investigation of Lagrangian dispersion is potentially useful for determining how clusters of passive or active particles e.g., oil, microplastic and marine species will spread from their initial positions (absolute dispersion) or from their center of mass (relative dispersion), and indeed how the cloud will disperse toward the zonal or the meridional directions, or otherwise showing how fast the cloud of particles will spread (relative diffusivity) depending on its initial release.
This paper presents, for the first time, the dispersion statistical properties and highlights the influence of the coherent structures in the occurrence of absolute and relative dispersion regimes in the BS using in-situ measurements and numerical model outputs.

Drifters
The BS drifter dataset used for this study was collected by 118 drifters deployed in 1999-2009. The drifters used are mainly composed of three different designs: (1) the Surface Velocity Program (SVP) drifters fixed at a depth of 15 m with a holey-sock drogue; (2) the Coastal Ocean Dynamic Experiment (CODE) drifters drogued at a depth of 1 m; and (3) the Compact Meteorological and Oceanographic drifters (CMOD), which include a 60 cm long aluminum cylindrical hull drogue with a sonobuoy case that is also fixed at 1 m depth [16]. The majority of collected data come from SVP drifters. The SVP drifters measure the near-surface currents in the first 15 m of the water column [2]. They are tracked by Argos satellite system, interpolated at 2-h intervals using the kriging method and finally filtered using a Hamming filter (cut-off period at 36 h) in order to remove higher frequency current components (tidal and inertial currents). Filtered positions were subsampled at 6-h intervals [2]. In this study, the drifters were handled in pairs (original and chance pairs) for the combination of the three drifter types (see for details, [23,24]). The pairs were selected for the initial separation distance D 0~3 .5-5 km, which is much smaller than R d (internal Rossby radius) ranging from 10 to 30 km for the BS [2,28,29]. The initial drifter pair tracks were considered every 5 days along the trajectories. This time interval is larger than the integral time scale estimated in the BS (T L~3 days) (see [2,23] for details).

Satellite Altimetry
The altimetry data used in this study is the gridded, delayed time products from AVISO (1999AVISO ( -2009, which is merged with the drifter data to estimate pseudo-Eulerian velocity statistics in bins of 0.25 • × 0.25 • × 1 day in order to present the mean BS currents and mean sub-basin eddies and gyres (for more details, see [16]).

Models
We built a high-resolution BS ocean model configuration with the 3D hydrodynamic Coastal and Regional Ocean COmmunity model (CROCO), which was developed around the Regional Oceanic Modeling System (ROMS) [30]. It resolves the hydrodynamic primitive equations and follows the Boussinesq and hydrostatic approximations, and has been designed to realistically resolve basin-scale, regional and coastal ocean processes at high-resolution [31]. The geographical domain of the configuration is extended from 26 to 42 • E and 40 to 47 • N (Figure 1), with a horizontal resolution of 1/32 • (~3.5 km).
The bathymetry was derived from the Global Earth Bathymetric Chart of the Oceans for the 2020 version (GEBCO_2020, Figure 1a) with a high spatial resolution of~0.45 km ( https://www.gebco.net (accessed on 10 October 2020)). Along the vertical water column, 20 vertical levels were parameterized using the σ-coordinate, with the thickness of the upper layer being around 2.3 m near the coast and 10.2 m in the central basin. In our configuration, the use of sigma coordinate allows providing a high resolution at the surface and the bottom layers in each point of the model grid. A similar study in the BS used the sigma coordinate scheme, which allows high resolution along the water column [32]. The surface layer of the simulation was forced using the National Center for Environmental Prediction (NCEP) atmospheric data with monthly averages and 1.9 • of spatial resolution [33,34]. As is done in [31,35], the atmospheric variables were used to provide the bulk conditions using a bulk formula representation of surface heat and freshwater fluxes [35]. The lateral Dardanelles boundary and initial oceanic conditions were set from monthly means of the Simple Ocean Data Assimilation (SODA) at 1/2 • of spatial resolution [36,37]. Internal tide can also play an important role in the dispersion near the coast and open ocean, especially, in shallow waters where the tide amplitudes become very important. In the case of the BS, the amplitude of the major tides is weak. Tidal sea level variations are only a few centimeters [38], compared to values larger than 50 cm in the Gulf of Gabes in the Mediterranean Sea [39]. In this configuration, river-runoff and tides are not parameterized.
The Lagrangian particle tracking model Ichthyop v3.3 is forced with velocity fields from the hydrodynamic model (CROCO) in order to simulate virtual particle tracks at intervals of 6 h. Ichthyop is a free Java tool designed to study the effects of physical and biological factors on the ichthyoplankton dynamics [40,41]. The Ichthyop model was recently used in the BS to simulate distribution and accumulation of floating marine debris [41]. Our model only considers the physical flow field generated by the CROCO model by omitting the biological parameters on ichthyoplankton dynamics [40][41][42]. In this study, we used the Lagrangian model to launch virtual drifter trajectories computed at each time step using a Runge-Kutta 4 forward scheme. We selected virtual drifter pairs from the virtual drifter trajectories for the same initial separation distance D 0~3 .5-5 km as the real drifter pairs displayed above by adapting the [24] approach.

Methods
The Lagrangian dispersion measures the spreading of particles as a function of time or separation distance. It is a powerful method to use both absolute (individual particles) and relative (particles pairs) dispersion to provide a statistical description of the turbulence statistics in the ocean.

The Turbulent Relative Dispersion Statistics
The relative dispersion characterizes the spreading of Lagrangian particles selected in pairs. It can be used as a fundamental method for the diagnosis of the patches amplification and mixing in the ocean.
To quantify the relative dispersion statistics, we defined the separation distance of the pair δ as the distance between two particles forming the pair with respect to their initial separation where r (1) and r (2) are the Lagrangian arc length of the particle trajectories forming the pair, t is the time and D 0 is the initial pair separation distance. The normalized fourth moment of the pair separation distances was defined as the kurtosis as a function of time [43] as follows: where < δ 2 (t, D 0 ) > is the relative dispersion and the average <.> is over all the drifter pairs with respect to their initial separation distance D 0 . Averaging the square of the pair velocity differences depicts the Second-order structure function where V (1) and V (2) are the Lagrangian velocities of two drifters separated with a predefined distance.
The relative diffusivity K(δ) which is the derivative of the relative dispersion, measures how fast the patches will be dispersed [44].
It can also be expressed with the Second-order structure function [45]: The time scale-dependent pair separation rate λ v (δ) was defined as follows [46]: The relative dispersion regimes are summarized in Figure 2a as suggested in [24,[43][44][45][46]. The well-known turbulent (a) relative and (b) absolute dispersion regimes, were t and T are the times, T L is the integral time scale, R d is the internal Rossby radius and δ is the separation distance between the pairs. The time parameter T allows estimating the theoretical exponential fit of Ku.

The Turbulent Absolute Dispersion Statistics
The absolute dispersion is defined hereafter as the measurement of the mean squared of the separation distance between individual particles as a function of time with particular attention to their initial positions. The absolute dispersion helps to measure the spreading of passive contaminant patches (e.g., heavy oil, pesticides, hydrocarbons, toxic metals) on a large scale in the ocean. The classical absolute dispersion approach answers the basic question of how the passive tracers spread in a turbulent field. It has been used in previous studies to distinguish the difference between the elliptic (strong rotation) and hyperbolic (dominant stretching) regions [23,24,47].
The absolute dispersion was defined as follow: where X(i, t) is the 2-D vector position of the drifter i at the time t, while X(i, t 0 ) is the vector position at t 0 . U is the mean velocities along the drifters tracks and <.> is the average over all individual drifters. At intermediate time scales, two anomalous absolute dispersion (5/3, elliptic) and (5/4, hyperbolic) power laws were previously found in a 2-D turbulent flow. These laws were detected from the calculation of the normalized Okubo-Weiss parameter (Q*) [47] when the particles move in elliptic (high vorticity gradient) and hyperbolic (high strain gradient) regions by using conditional averages. This indication implies that the anomalous regimes are related to the flow topology (hyperbolic and elliptic areas). where: The elliptic regime is related to the regions where the rotation is stronger than the deformation ( s 2 ∼ 0 , then Q * ~− 1), whereas the hyperbolic regime is related to domains characterized by dominant strain gradients ( w 2 ∼ 0, then Q * ~1 ). There are also regions where the strain and the vorticity are balanced (s 2 ∼ w 2 , then Q * ∼ 0) [47][48][49][50]. The absolute dispersion regimes are summarized in Figure 2b, as suggested in [23,24,47,50].

Circulation of the BS
The mean surface circulation in the BS was retrieved using the method of [16]. Drifter measurements and satellite altimetry data sets were combined for the period 1999-2009 to estimate the mean geostrophic circulation field and to describe the surface sub-basin scale coherent structures and currents ( Figure 1b). The surface circulation is strongly characterized by the RC with the presence of gyres and eddies ( Figure 1b). The winddriven RC is very persistent; its high meandering signature close to the coasts leads to the generation of sub-basin scale and mesoscale structures such as the EG, WG, as well the edges of the BE and SevE [1,2,16]. According to the results of [16], the mean currents were generally stronger in winter, mainly in the years 2002-2006, and the most intense activity of the BE occurred in summer in 2006. We selected two months, January and August that are representative of winter and summer conditions, in two different years (2002 and 2006). We performed a qualitative and quantitative comparison between current fields obtained by merging drifters and altimetry (Figure 3a  The seasonal goodness of fit of the model can be represented by the RMSE = (ud − um) 2 + (vd − vm) 2 , where ud and um are the zonal velocity compo-nents for the drifter-altimetry and model, respectively, while vd and vm are their meridional components; the discrepancy between the drifter-altimetry and the model derived geostrophic velocities was found to be 3.67 cm/s and 3.24 cm/s in winter, whereas in summer, it was found to be 4.38 cm/s and 3.83 cm/s in 2002 and 2006, respectively. The discrepancy between observed and modeled currents is larger in summer, in particular in summer 2002, probably due to the low-resolution of atmospheric forcing (NCEP) used for the simulation during this season, and/or from the reduced freshwater flux from river inputs.
In winter, the RC is very persistent and tends to form a single cyclonic loop along the coasts, whereas in summer the flow exhibits a meandering signature.     Larger values of KE (of~0.03 m 2 s −2 ) were observed in winter in the branches of the RC located along the Anatolia coast, the Crimea peninsula, and in the western corner of the BS (Figures 4a and 6a), in agreement with previous studies [2,16]. In winter, the strength of the RC in the intermediate layer is reduced with respect to the surface, and the signature of SevE and BE eddies are clearly identified (Figures 4b  and 6b). Another anticyclonic mesoscale eddy, known as SakE [16], is located along the south-western coast and centered between 30 • E and 31 • E. In the deep layer, the signature of the RC disappears, and the model derived currents show a couple of coastal anticyclonic eddies (Figures 4c and 6c).
With a reduced cyclonic basin-scale circulation in the interior of the BS, it is possible to discern the signature of numerous forms of mesoscale eddies circulation located along the coast. In summer, the signature of the RC is less evident in the intermediate layer; in particular in August 2002, the main signature was limited to the EG (Figure 5b); numerous anticyclonic mesoscale structures were retrieved along the coast (Figures 5b and 7b).
In the deep layer, the signature of the RC is not apparent and the coastal circulation is characterized by a weak anticyclonic flow, with larger speeds observed in the regions of mesoscale coastal eddies.
The normalized Okubo-Weiss parameter Q* Equation (8) computed from the model is shown in Figures 4-7 (right panels) for similar scenarios as shown in Figures 4-7 (left panels) discussed above. For Q* close to 1 (hyperbolic regions) or −1 (elliptic regions) two tracers move away from each other, whereas for Q* equivalent to zero, the tracers stay where they are located. Hyperbolic regions (stretching, red color) were detected outside the eddy, whereas elliptic regions (vorticity, blue color) were detected in their inner part. The spatial distribution of Q* parameter derived by the model (the right panels of Figures 4-7) shows the presence of numerous mesoscale structures and filaments, quite similar to the circulation shown by the black arrows highlighted (in the relative left panels of Figures 4-7). In winter, the surface layer is dominated by hyperbolic regions, related to the strengthening of the RC (Figures 4d and 6d). In the intermediate layer, the Q* clearly describes two large elliptic regions, corresponding to the EG and WG; the area between these two sub-basin gyres and the coast is characterized by numerous mesoscale structures (Figures 4e and 6e). Mesoscale eddies in the BS have a life time larger than four weeks and diameters larger than 40 km [52] and are therefore well illustrated by the monthly averages of Figures 4-7. In the deep layer, the signature of sub-basin gyres is less evident and the entire basin is marked by smaller scale structures (Figures 4f and 7f). In the surface layer, the elliptic regions are more commonly detected in summer than in winter, due to the meandering flow of the RC and to the marked presence of mesoscale eddies (Figures 5d and 7d). The signature of these mesoscale eddies also extends into the intermediate layer (Figures 5e and 7e). In the deep layer, mesoscale eddies persist in the central area, whereas numerous mesoscale structures are observed in the coastal regions (Figures 5f and 7f).

Dispersion
The mesoscale features described in the Section 3.1 influence the water mass transport and, in turn, the biodiversity of the BS. In this context, the description of the dispersion processes associated with these coherent structures becomes very interesting. Lagrangian dispersion estimated using drifter tracks and/or model are presented for the surface, intermediate, and deep layers. The effects of the RC and eddy structures on the variability of dispersion in the water column are discussed.

Relative Dispersion Surface Relative Dispersion
The larger part of the Lagrangian drifter trajectories deployed in winter/spring ( Figure 8a) follows a single cyclonic loop which is the signature of the RC, whereas in summer/fall, they fluctuate more due to the presence of meandering currents (Figure 8d). Drifter tracks display the strong signature of the cyclonic and anticyclonic eddies already reported in a previous study [2]. The results of the Lagrangian particle tracking model for the January and the August of 2002 and 2006 are shown in Figure 8b,c and Figure 8e,f, respectively. The 30-days virtual drifter trajectories show that in winter (January) the drifters tend to form a single cyclonic loop which is the edge of the RC, while in summer (August) the flow becomes more meandering in accordance with the results shown by drifter tracks (Figure 8a,d). In order to show the presence of the statistical relative dispersion regimes, we discuss the drifter pair dispersion statistics by comparing three different accurate metrics, namely the relative diffusivity (Figure 9; left panel), the timescale dependent pair separation rate as a function of pair separation distances (Figure 9; right panel), and the fourth moment of pair dispersion distances (Kurtosis) as a function of time (Figure 9; insets). These metrics are derived from the Lagrangian pairs, as shown in Figure 8, for both drifters and model data. For scales of 10 to 100 km, the relative diffusivity K(δ) and the time-scale λ v (δ) curves show the presence of the Richardson regime in all cases because the pair separation distances are comparable to the scales of dominant eddies, in agreement with the results of [23,24]. The Richardson regime occurred for spatial scales between 10 and 100 km due to the signature of an inverse energy cascade range. This result is supported by the Kurtosis (Ku), where Ku is about 5.6 for intermediate time scales of 7-15 days. In contrast, at large scales (100-300 km), the pair statistics show the occurrence of the diffusive-like regime only for the case of drifter pairs selected in winter/spring, when Ku is close to 2 for time scales larger than 10 days. In fact, the occurrence of the diffusive regime in this scenario implies that the majority of drifter pair velocities are uncorrelated, and the pairs lose memory of their initial separation distances, as suggested in previous studies [23,53]. For separation scales below 10 km, the model only underlines the occurrence of the non-local exponential regime, related to an enstrophy cascade range, in August 2002, while in the other cases, this regime is absent. These numerical results are in agreement with drifter derived statistics.

Intermediate and Deep Relative Dispersion
Until now, only a few studies have focused on the dispersion in the intermediate and deep layers of the ocean [53][54][55][56]. In this study, a powerful tool such as a high-resolution numerical model can be employed to investigate the Lagrangian dispersion along the water column. In this context, we highlight how the deep currents and coherent structures could influence the dispersion. Similar to in Figure 8 The virtual drifter pairs were obtained from the same initial virtual drifters that were launched in the surface layer to compare the drifter trajectory behaviors (for details, see Figure 8 vs. Figure 10). Relative dispersion statistics, estimated using virtual drifter pair trajectories in the intermediate and deep layers, are shown in Figure 10 and displayed in Figure 11 (2002) and Figure 12     The difference in the occurrence of the non-local regime is related to the different dynamics and eddy properties in each layer of the water column and to the seasonal variability that induces the pairs to escape or be trapped in fronts characterized by a high KE gradient (see Figures 4-7). The exponential fits presented in the insets of Figure 9, Figure 11, and Figure 12 were obtained from the following theoretical expression Ku = exp(8t/T), as shown in (Figure 2a). For the exponential model, i.e., at fixed time te = 3 days, T = 8te/ log[Ku(te)] (in days); for more details, see [24,43].

Absolute Dispersion Surface Absolute Dispersion
The surface absolute dispersion was analyzed using drifters in winter/spring and summer/fall; it was also modeled for January/August of 2002 and 2006, as shown in Figures 13 and 14.  Numerous drifter pairs were selected in a spatial range between 0 and 100 km to guarantee a consistent number of segments (individual drifters) for the absolute dispersion calculation (insets of Figure 13) and thus ensure a high statistical robustness (see [24] for more details), while the absolute dispersion for the model was estimated following the same drifter pairs that were used in the relative dispersion statistics (Figure 8) due to the short life times of the virtual drifters (30-days drifter trajectories). The absolute dispersion was found to be anisotropic, with the zonal component being much larger than the meridional one in all cases due to the strong effect of the RC flowing mostly in the zonal direction. At small (<3 days) and large (>10 days) time scales, quasi-turbulent absolute dispersion regimes ballistic and random-walk, respectively, were detected in all cases. Examining the surface absolute dispersion behaviors at intermediate times, this reflects the occurrence of the two anomalous regimes: elliptic (5/3) and hyperbolic (5/4). The hyperbolic regime was easily detected in winter/spring and January 2002 and 2006 (Figures 13a and 14a,c), while the elliptic regime was found in summer/fall and August 2002 and 2006 (Figures 13b and 14b,d) using both drifters and the model.
The two anomalous regimes were supported by showing the temporal evolution of the absolute dispersion curves normalized by t 5/4 (Figures 13c and 14e,g) and t 5/3 (Figures 13d and 14f,h), for the hyperbolic and elliptic regimes, respectively. In fact, a quasi-plateau was found, reflecting the dominance of the hyperbolic and elliptic regimes in any cases discussed above. The occurrence of the hyperbolic regime over a relatively cold season (winter/spring and January) is related to the dominance of hyperbolic regions, which were detected by the Q* parameter; more than 70% of hyperbolic areas were detected in this scenario (see the insets of Figures 4d and 6d). The (5/3) law was only observed in relatively warm seasons (summer/fall, August), essentially due to the fact that particles were launched into the elliptic regions for a persistent meandering flow (blue regions, Figures 5d and 7d).  Figure 16f,h for the hyperbolic and elliptic regimes, respectively. The quasi-classical absolute dispersion regimes (ballistic~t 2 and random-walk~t) were also detected at depth.

Discussion
A general description of the Lagrangian absolute and relative dispersion statisticswas carried out using Lagrangian drifter tracks and a high-resolution numerical model for the time period of 1999 to 2009. The numerical simulation results presented in this research are in agreement with previous studies [1,2,16,[57][58][59][60][61] and provide an overview of the mean sub-basin and mesoscale features of the surface, intermediate, and deep layers of the water column. The model resolves the surface circulation and the mean mesoscale features in the two years selected (2002 and 2006) during the coldest (January) and the warmest (August) months. The RC is very persistent, especially in January in the surface layer, and is very strong near the coast, while at depth it becomes less pronounced due to the absence of the wind stress influence. In summer (August), the RC signature is relatively weaker and the flow exhibits a meandering nature, probably due to mixing and baroclinic instability. An underestimation of the simulated surface circulation was observed in both seasons, with larger discrepancies in summer (Figure 3f,h), compared to the flow field shown by drifters and altimetry (Figure 3b,d). This discrepancy could possibly be related to the low-resolution of the atmospheric NCEP dataset used to force the model and to the weaker freshwater flux observed in summer (in particular in August 2002; not shown).
We used the normalized, model derived Okubo-Weiss parameter Q* to show that the BS mesoscale eddies can be classified as elliptic (near to the eddy cores) and hyperbolic (eddies surrounding) regions. Fluctuating from −1 to 1, the Q* detects the elliptic regions when it is close to −1, while the hyperbolic regions were found for Q*~1. Close to zero, the elliptic and hyperbolic regions are balanced and the flow is more homogeneous; this condition was well detected at depth.
The relative dispersion plays an important role on the turbulence statistics, because it gives a general idea of water mass spreading in the ocean and answers to complex phenomena such as chaotic turbulent flow and the contribution of currents and eddies to the transport process. The effect of mesoscale structures on the relative dispersion regime is well described from the occurrence of the surface Richardson regime (Figure 2a) derived from drifters and the model simulation. The Richardson regime takes place at intermediate space scales. It was observed for a separation distance δ~1-3 km in the Gulf Stream [20] and for scales smaller than 10 km in the northern Gulf of Mexico [46], while it extends to 150 km in its southern part [43]. In the Mediterranean sub-basins, it occurs when the pair separation distances become comparable to the scales of the main eddies [24]. In the BS, this regime was detected for δ ∼ R d due to an inverse energy cascade range. Previous studies demonstrated that the presence of the shear/ballistic regime is essentially due to the strong boundary currents and fronts representing a high KE [62]. This intermediate regime was found in previous studies e.g., recently in the Mediterranean Sea, it was only detected in the Adriatic sub-basin due to the presence of the coastal Adriatic jets close to the large boundaries (for more details, see [24]); this condition is absent in the BS. However, this regime is absent because the pairs would escape from fronts characterized by a high KE and move toward regions with less energy gradients.
In general, the non-local exponential regime was detected for spatial scales smaller or comparable to R d , from model results (Figure 2a). The presence of this regime can be related to the effect of coherent structures on the pair separations at small spatio-temporal scales rather as the ability of the model to solve the exponential regime for small scales. This is because this regime is related to the role of the coherent structures at scales less than R d . According to these results, it is clear that the model is able to resolve the complex turbulent variability and could also help to explore the non-local regime along the water column (Figures 11 and 12).
The relative diffusivity in winter/spring/January is larger than in summer/fall/August (Figures 9, 11 and 12), due to the influence of the strong RC. In this condition, the passive patches (for e.g., pollutants) will disperse much faster in winter/spring/January. In summer/fall/August, the weak diffusivity is related to the meandering flow, the eddies trap drifters over a few weeks, and the fact that the pairs are not dispersed over large distances.
The exponential regime, which was detected only by the model for scales <10 km, has been an open question in previous works [19][20][21][22][23][24]. This regime is hardly detected by in-situ measurements at small spatial scales, due to the non-uniform nature and no homogeneity of the drifter data. Drifters were released indifferent geographical areas characterized by different dynamics. As a result, this nonlocal regime depends strongly on the initial release. In our study, the presence of this regime is detected due to the relative diffusivity for intermediate space scales, and is well supported by Ku and λ v (δ) as a function of time and spaces, respectively. On the other hand, the Richardson regime can extend to small scales, indeed this does occur with the relative diffusivity increasing quickly and inhibiting the occurrence of the exponential regimes. A similar scenario for the faster increase of diffusivity inhibiting the presence of the non-local exponential regime is shown in [23]; this condition is absent in our statistics.
The anomalous absolute dispersion regimes, hyperbolic (5/4) and elliptic (5/3) (Figure 2b) regimes were detected in the surface layer and as a function of depth in all cases (Figures 13 and 14). The hyperbolic regime is related to the topology of the BS. It was detected only in winter/spring (Figures 13-16). Meanwhile, the elliptic regime, related to the meandering pathway of the RC and to the mesoscale activity, was detected in summer/fall. Results derived from the estimation of Q* (Figures 4-7) show predominant hyperbolic areas (Q* > 0) in the surface layer, both in winter (more than 70%; Figures 4d and 6d), in good agreement with the absolute dispersion results (Figures 13 and 14), and in summer (~60%; Figures 5d and 7d), while the absolute dispersion shows a prevalent elliptic regime (Figures 13 and 14). The two results do not contradict each other, as the estimate of Q* percentage refers to the whole basin (where the hyperbolic regime is clearly dominant), whereas the absolute dispersion refers to the specific paths followed by the real and simulated drifters, most of which are trapped within the eddy inner parts in summer, emphasizing the influence of the elliptic areas (Figure 8d-f). A similar scenario was also observed in the intermediate and deep layers.

Conclusions
This paper provides, for the first time, a general overview of the role of the persistent currents and eddies in the BS for the occurrence of both relative and absolute dispersion regimes. In terms of the relative dispersion (for intermediate time and space scales), at the surface layer (15 m), the flow implies the occurrence of the Richardson regime due to the presence of an inverse energy cascade range. For short times and small separation distances, the relative dispersion statistics denote the existence of the non-local exponential regime where the pair separations are strongly driven by eddies when their scales are larger than the pair separation distances. In the intermediate (150 m) and deep (750 m) layers, the exponential and Richardson regimes occurred due to the effect of the intermediate and deep BS circulation features.
In terms of absolute dispersion regimes, the results investigate the important role of the BS topology on the occurrence of the hyperbolic regime in winter/spring and January, and on the occurrence of the elliptic regime in summer/fall and August. The normalized Okubo-Weiss parameter shows for the first time the spatial distribution of the elliptic and hyperbolic regions in the BS for high vorticity gradients and dominant stretching, respectively. The absolute dispersion of BS is anisotropic, in agreement with the result of [2]. This paper outlines how a turbulent flow field could induce the discrepancy between the theory, the reality, and the simulation in one of the most isolated areas in the global ocean. A high-resolution ocean model could be a powerful tool to further explore this discrepancy. However, the model simulation is always limited by the resolution. For a more realistic understanding/predictions, it would be interesting to develop a simulation with more accurate boundary and surface conditions, which seems to be critical for the representation of the flow field.