Influence of the summer deep-sea circulations on passive drifts among the submarine canyons in the northwestern Mediterranean Sea

Marine biophysical models can be used to explore the displacement of individuals in and between submarine canyons. Mostly, the studies focus on the shallow hydrodynamics in or around a single canyon. In the northwestern Mediterranean Sea, knowledge of the deep-sea circulation and its spatial variability in three contiguous submarine canyons is limited. We used a Lagrangian framework with three-dimensional velocity fields from two versions of the Regional Ocean Modeling System (ROMS) to study the deep-bottom connectivity between submarine canyons and to compare their influence on the particle transport. From a biological point of view, the particles represented eggs and larvae spawned by the deep-sea commercial shrimp Aristeus antennatus along the continental slope in summer. The passive particles mainly followed a southwest drift along the continental slope and drifted less than 200 km considering a pelagic larval duration (PLD) of 31 d. Two of the submarine canyons were connected by more than 27 % of particles if they were released at sea bottom depths above 600 m. The vertical advection of particles depended on the depth where particles were released and the circulation influenced by the morphology of each submarine canyon. Therefore, the impact of contiguous submarine canyons on particle transport should be studied on a case-by-case basis and not be generalized. Because the flows were strongly influenced by the bottom topography, the hydrodynamic model with finer bathymetric resolution data, a less smoothed bottom topography, and finer sigma-layer resolution near the bottom should give more accurate simulations of near-bottom passive drift. Those results propose that the physical model parameterization and discretization have to be considered for improving connectivity studies of deep-sea species.

1 Introduction

Lagrangian particles coupled with three-dimensional hydrodynamic models are useful to assess the impact of ocean circulation on the drift of small elements or individuals. It allows for the exploration of various scenarios of dispersal and increases knowledge in several fields of marine systems, such as predicting the direction of oil spills (Jones et al., 2016); understanding the circulation of microplastics (Lebreton et al., 2012); estimating the impact of aquaculture (Cromey and Black, 2005); and describing the spatial dispersion of crustaceans, fish, and other marine larvae (Ahumada-Sempoal et al., 2015; North, 2008; Ospina-Aìlvarez et al., 2013). Particles in the open ocean are susceptible to advection between locations, influenced by regional currents and mesoscale features such as eddies and meanders (Ahumada-Sempoal et al., 2013). The particles without implemented behavior are called passive and are efficiently used to study the inactive transport of elements. They also provide a useful approximation of larval drifts when the ecological knowledge on the early life cycle is scarce. The individual-based models (IBM) configure each Lagrangian particle with characteristics parameterized by the modeler. In marine ecological studies, IBMs provide a representation of the potential connectivity between geographically separated subpopulations, with implications for fishery management and conservation plans (Andrello et al., 2013; Basterretxea et al., 2012; Kough et al., 2013; Siegel et al., 2003).

Simulation of marine passive drifts provides a picture of the highly dynamic circulation in the vicinity of the submarine canyons. Influence of canyons on the transport of particles and living organisms has begun to raise interest, but only a few studies are dealing with the Lagrangian transport of particles specifically in those small-scale topographic structures (Ahumada-Sempoal et al., 2015; Kool et al., 2015; Shan et al., 2014). In the northwestern (NW) Mediterranean Sea, several submarine canyons, whose heads are incising the continental shelf (i.e., Cap de Creus, Palamós, and Blanes canyons), generate mesoscale features from the main circulation pattern (Millot, 1999) called the Northern Current (NC). When the NC crosses a canyon, it causes a downwelling on the upstream wall and an upwelling on the downstream wall in the 100–200 m layer depth (Ahumada-Sempoal et al., 2015; Jordi et al., 2006), or it causes the opposite at lower depths (Flexas et al., 2008). Within the canyons, near-bottom currents produce a closed circulation (Palanques et al., 2005; Solé et al., 2016), which loses strength with the depth (Flexas et al., 2008; Granata et al., 1999). The canyon shape enhances the downwelling of water and its components (sediment and organic carbon; Puig et al., 2003, biogenic or inorganic particles; Granata et al., 1999), and it provides favorable and diverse habitats for benthic species (Fernandez-Arcaya et al., 2017). Those exchanges proceed at different velocity rates following the condition of the waters. For example, cascading of winter water masses drives suspended substances along the funnel structure of the canyons (Canals et al., 2009). In summer, the stratification of the water column is well-established in the NW Mediterranean Sea decoupling the mixed layer from the rest of the water column (Rojas et al., 1995) and slowing the downward sink of biogenic particles (Granata et al., 1999). During that season, the circulation of the NC is shallower (250 m deep), wider (50 km wide) and less intense than in winter (Millot, 1999).

To date, the impact of the circulation on the dispersal of particles was estimated with hydrophysical data (Ahumada-Sempoal et al., 2015; Granata et al., 1999; Rojas et al., 1995) and with sediment traps (Durrieu de Madron et al., 1999; Lopez-Fernandez et al., 2013) in shallow waters above the NW Mediterranean submarine canyons. None of the previous studies simulated the particle transport at a deeper layer or showed the particle exchanges among different water layers of submarine canyons. Information on the roles of the upper part of the canyon on both horizontal and vertical transport of particles is not clear (e.g., 0–300 m; Granata et al., 1999), and for deeper parts this information is rarely provided. In accordance with those facts, we questioned if the deep-sea circulation mechanism in each submarine canyons of the NW Mediterranean Sea can be generalized.

In this study, deep passive drifts along the NW Mediterranean continental slope are simulated to determine the influence of the summer circulation modulated by the presence of the submarine canyons. We partly focused on the drift sensitivity to deep circulation using the outputs of two hydrodynamic numerical models to approach the transport uncertainties. With this study, we are expecting to provide a first insight into the deep connectivity between submarine canyons to improve the management of their grounds, namely of the deep-sea shrimp Aristeus antennatus.

2 Material and methods

The influence of deep circulation on passive drifts among the submarine canyons of the NW Mediterranean Sea was analyzed through the dispersal of Lagrangian particles.

2.1 Hydrodynamic models

The sensitivity of dispersal to the circulation was generated with the use of two available hydrodynamic models, which covered the area of the submarine canyons in the northwestern Mediterranean Sea (see Fig. 1). Because of the parameterization related to the topography (i.e., bathymetric data sets and smoothing) and the vertical layer discretization in the two hydrodynamic models, small differences in morphology of the submarine canyons were expected to impact the flow circulation.

Figure 1 Map of the study area. Release zones are represented along the northwestern Mediterranean Sea: open slope (OS) 1, Cap de Creus canyon (CCC), OS 2, Palamós canyon (PC), OS 3, Blanes canyon (BC), and OS 4. The red arrow indicates the direction of the Northern Current. Bathymetry was extracted from EMODnet (http://www.emodnet-bathymetry.eu, last access: December 2018).

The two hydrodynamic models were climatological simulations using the Regional Ocean Model System (ROMS; Shchepetkin and McWilliams, 2005), a free-surface terrain-following primitive equations ocean model. ROMS includes accurate and efficient physical and numerical algorithms (Shchepetkin, 2003). The two models were referred to by their implementation versions, ROMS-Rutgers and ROMS-Agrif (Debreu et al., 2012), for the sake of clarity. In the present paper, the ROMS-Rutgers configuration and its validations are presented, while ROMS-Agrif has already been used and validated in Ahumada-Sempoal et al. (2013, 2015).

The ROMS-Rutgers model was forced by a climatological atmospheric forcing. The air temperature, shortwave radiation, longwave radiation, precipitation, cloud cover, and freshwater flux used to force the model came from the ERA-40 reanalysis (Uppala et al., 2005). Surface pressure came from the ERA-Interim reanalysis (Dee et al., 2011). All these variables had a spatial resolution of 1^∘ and a time resolution of 6 h. QuikSCAT blend data were used for wind forcing (both zonal and meridional). The wind had a spatial resolution of 0.25^∘ and a time resolution of 6 h. The boundary conditions were obtained from NEMO (available from https://www.nemo-ocean.eu/, last access: March 2012; these simulations were reported in Adani et al., 2011) and interpolated to the ROMS grid to define a sponge layer of 10 horizontal grid points with a nudging relaxation time of 30 d. This methodology for the implementation of the model in the area followed the same procedure as the one already tested in the Alboran Sea (Solé et al., 2016). This model implementation was already used in previous publications as a hydrodynamic model coupled with a fisheries model (Coll et al., 2016).

The simulation domain ranges from 0.65^∘ W to 6.08^∘ E and from 38 to 43.69^∘ N. The grid resolution is 2 km (with 256×384 grid points horizontally) and the vertical domain is discretized using 40 vertical levels with a finer resolution near the surface (S coordinate surface and bottom parameter are θ_s=5 and θ_b=0.4). Thus, the thickness of the near-bottom layer on the continental slopes delimited by the 500 and 800 m isobaths is 24 m (±3.2 m). The advection scheme used in our simulations was MPDATA (recursive flux-corrected 3-D advection of particles; Smolarkiewicz and Margolin, 1998) and Large–McWilliams–Doney (LMD) mixing as a subgrid-scale turbulent mixing closure scheme (Large et al., 1994), also known as the K-Profile Parameterization (KPP) scheme. The air–sea interaction used for the boundary layer in ROMS was based on the bulk parameterization (Fairall et al., 1996). Bathymetry grid with a horizontal resolution of 100 m was provided by the Geoscience department of the Institute of Marine Sciences (ICM-CSIC, Barcelona) and fitted to the grid of the model. We ran the ROMS-Rutgers model version 3.6 using climatological atmospheric forcing and boundary conditions. The initial conditions to start the simulations were obtained using the same interpolated fields as the ones used for the boundary conditions for all variables. After an 8-year spin-up period with a baroclinic time step of 120 s, we used the ninth year as the study period.

The regional ROMS-Rutgers implementation has been daily saved and validated. First, we compared the simulations with the coarser model used for the initial and boundary conditions (NEMO). Second, the model was compared to the ARGO float vertical data profiles of temperature and salinity to test the correct structure of the simulated water column within the year. For the ARGO tests, we selected 1900 casts in the area from October 2003 to December 2012. The data collected cover depths from the surface down to a maximum of 720 m. The ARGO vertical profile resolution was 5 m from the surface to 200 and 25 m from 200 to 720 m. We grouped the ARGO casts by months and by six subareas dividing the domain according to the Coriolis force divided by the bottom depth (Fig. 2). Then, we selected the subareas that have more than 30 ARGO profiles for each month, we calculated a monthly averaged profile of ARGO, and we compared them with the modeled climatologic profile (monthly average).

Figure 2 Map of the ROMS-Rutgers model domain and position of ARGO buoys. The color bar represents the Coriolis parameter (f) divided by the depth (h) (f∕h in a zebra color palette). The six subareas selected to validate the model using ARGO float observations are shown by yellow polygons. The black dots in the domain represent all the ARGO profiles considered.

The main general behavior of the ROMS-Rutgers model simulation was coherent with the lower-resolution model (NEMO) and with the reported hydrography of the northwestern Mediterranean Sea. The model also successfully reproduced the main seasonal behavior of the different water masses. A Taylor diagram (Taylor, 2001) was used to display the correlation, root mean square error (RMSE), and standard deviation (SD) between the monthly averages of the ROMS-Rutgers model profiles and ARGO profiles by subareas (Fig. 3). The comparison showed reasonably good correlations of statistical significance. For temperature, the correlation between model and ARGO was higher than 0.7 during 10 of the 12 months, while for salinity it was higher than 0.95 for all the months. As it is shown in Fig. 3, during the month of September, the correlations in the Taylor diagram for temperature and salinity were higher than 0.9 in both variables. The month of September was selected because the thermocline in our domain disappears and the mixing process at that period of the year is likely to vary the most. Consequently, it shows that the circulation over different depths has been well reproduced in the model.

Figure 3 Taylor diagrams for validation of ROMS-Rutgers. Example of comparison between buoys (ARGO floats) and ROMS-Rutgers model simulations during September. (a) The temperature of float versus model and (b) the salinity of float versus model. Radial from the origin (0, 0) are the correlation coefficients, internal circles centered on the origin are standard deviations, and internal circles centered at the buoy position are root mean square differences (RMSDs).

The second set of velocity and thermodynamic fields was provided by ROMS-Agrif, built and validated in Ahumada-Sempoal et al. (2013). The simulation domain ranges from 40.21 to 43.93^∘ N and 0.03^∘ W to 6^∘ E. It has a finer horizontal resolution (∼1.2 km), with 32 sigma levels of a finer resolution near the surface and coarser resolution near the bottom (S coordinate surface and bottom parameters are θs=7 and θb=0) than ROMS-Rutgers. The average thickness of the near-bottom layer is 54 m and is approximately 2 times thicker than the near-bottom layer from ROMS-Rutgers. The model was one-way nested from a coarse regional resolution model of 4 km. Bathymetry was derived from the ETOPO2-arcminute model (with horizontal grid resolution of 3 km) from Smith and Sandell (1997). Surface forcing and initial and boundary conditions were built with the ROMSTOOLS package (Penven et al., 2008), running a 10-year simulation of the model.

Figure 4 Vertical velocity and horizontal velocity vectors at bottom in the northwestern Mediterranean Sea in the (a) ROMS-Rutgers and (b) ROMS-Agrif models. Colors from blue to red according to the downward or upward (negative or positive) vertical velocity. Continuous black lines are the isobaths adapted from the bathymetry of the hydrodynamic models.

In the ROMS-Rutgers configuration, the bottom circulation followed a southward direction along the bottom slope, except in the area south of Blanes canyon (Fig. 4). Over the bottom floor of the continental slope, the highest velocities reached 8.5 cm s⁻¹ on the southern mouth of Palamós canyon. The fastest vertical velocities were at sea bottom deeper than 400 m and are mostly located among the submarine canyons. The minimum and maximum velocity values (−1.4 and 1.7 mm s⁻¹, respectively) were located in Palamós canyon. In ROMS-Agrif, the hydrodynamic model had higher near-bottom currents in the areas at the north and south of Palamós canyon. The maximum velocity reached 11.6 cm s⁻¹ off the continental slope between Palamós and Blanes canyons. Higher intensities of vertical near-bottom current were also in those areas (see Fig. 4) even though the most extreme values of vertical velocity were −1.23 and 1.11 mm s⁻¹ on the bottom off the continental rise.

2.2 Practical study

In this study, we considered some characteristics of the deep-sea shrimp Aristeus antennatus for the initial configuration of the simulated passive drifters. A. antennatus is a benthic commercial species distributed over the area of study (Demestre and Lleonart, 1993; Sardà et al., 1997), and it is an important source of income for the fishing harbors of the NW Mediterranean Sea (Gorelli et al., 2014).

The spawning peak of A. antennatus – and release of particles – occurred in summer (July to September; Demestre and Fortuño, 1992) in the sea bottom delimited by the isobaths of 500 and 800 m, i.e., on average 650 m (Sardà et al., 1997; Tudela et al., 2003) at midnight (Schram et al., 2010). The ecology of the early pelagic stages (eggs and larvae) is hardly known due to the difficulty to catch the larvae in the open sea and to keep the adults in captivity alive.

For setting the duration of simulated drift, we used the pelagic larval duration (PLD) of the shrimp, i.e., the time a larva spends in the water column from spawning to the first post-larval stage. As individuals of the superfamily Penaeidae (which contains the species A. antennatus) release eggs in the environment before hatching (Dall, 1990), the PLD definition was extended to account for the duration of the embryonic stage. To overcome the situation of an unknown PLD for A. antennatus, we fitted a linear model to the relation between the duration of eggs developing into the post-larvae stage from 43 penaeid species by reviewing research articles (see Table S1 in the Supplement) and the temperature of the water in which the larvae were reared. The linear model, whose initial assumptions were verified (see Fig. S1) and had a coefficient of determination R²=56 %, was the following:

$$\begin{array}{}\text{(1)}& \mathrm{PLD}=\mathrm{64.71}\cdot \exp (-\mathrm{0.06}\cdot T),\end{array}$$

where T is water temperature in degrees Celsius. Then, the effect of water temperature on the duration of the larval stage had the shape of an exponential law (O'Connor et al., 2007; Pepin, 1991). When particles were at the deep-sea bottom, we could estimate their PLD to be 31 d because the seawater temperature was approximately 13.2 ^∘C. The duration appeared to be shorter than the PLD for other deep-sea species (Arellano et al., 2014: Etter and Bower, 2015; Young et al., 2012) but was in agreement with PLD for temperate invertebrates (Levin and Bridges, 1995; Thatje et al., 2005; Williamson, 1982) and PLDs of penaeid species predicted in Dall (1990).

Using the passive drift simulations, we started to explore the simplest hypothesis around the larval ecology, which corresponds to the absence of larval behavior (such as egg buoyancy or swimming abilities). The drift of particles is used as a first step to approach the unknown larval cycle. The good management of the fishery has raised the interest of science to have a better knowledge of the species ecology. Last, the particles are used to estimate potential connections among the fishing grounds of A. antennatus. It contributes to complement the research made in parental genetic studies (Sardà and Company, 2012), which have shown the high dispersive capacity of A. antennatus without revealing its connectivity paths.

2.3 Biophysical model

Particle dispersal was simulated with the open-source IBM code Ichthyop (Lett et al., 2008) version 3.3, which is a free user-friendly toolset used in numerous individual dispersal studies in the western Mediterranean Sea (e.g., Ospina-Aìlvarez et al., 2013; Palmas et al., 2017). Displacements of virtual particles are computed by integrating the differential equation using a Runge–Kutta 4th order (RK4) advection scheme:

$$\begin{array}{}\text{(2)}& \mathrm{d}\mathit{X}=\mathit{U}\cdot \mathrm{d}t,\end{array}$$

where dX is the three-dimensional displacement vector of the particles during a time step dt of 30 min under the velocity vector U from the hydrodynamic models. RK4 is a stable and reliable multistep method for numerical integration (North et al., 2009; Qiu et al., 2011) and is commonly used in Lagrangian dispersal models (North et al., 2009; Chen et al., 2003; Paris et al., 2013). The Ichthyop algorithm used trilinear and linear interpolations in space and time from the daily average ROMS velocity field output to each particle position at all time steps. Attributes (geographical coordinates and depth) of each particle were saved in an output file at a daily time step.

To reach the optimal number of released particles, we carried out an analysis following a slightly modified protocol proposed in Hilário et al. (2015) and Simons et al. (2013). The drifts of 200 000 advective-only particles were carried out over 31 d (PLD computed in the above section), released at midnight for each day of July and August at the bottom. A subset of N particles (N=1000, 2000, 5000, 10 000, 25 000, 50 000, 75 000, 100 000, and 150 000) was randomly chosen among the 200 000 originally released. Then, relative to the N particles released, the density of particles in each cell of a two-dimensional (2-D) horizontal grid with a 2 km resolution (i.e., density matrix) was computed.

This operation was replicated 100 times, leading to a sample of one hundred 2-D horizontal density matrices for a given N. The difference between the 100 replicates of that sample was evaluated by calculating the two-by-two fraction of unexplained variance (FUV):

$$\begin{array}{}\text{(3)}& \mathrm{FUV}=\mathrm{1}-r^{\mathrm{2}},\end{array}$$

where r² is the squared Pearson linear correlation coefficient within the density matrices. Based on the results, the number of particles showing a FUV lower than 5 % rounded at 20 000 units. This number was set as the reference number of particles to be released per event.

2.4 Dispersal analysis

Seven zones were defined along the NW Mediterranean Sea continental slope including the Cap de Creus, Palamós, and Blanes submarine canyons (CCC, PC, and BC, respectively) and the open slopes (OS) between them (see Fig. 1). Those zones were also drawn to incorporate the continental shelf, slope, and rise as large polygons. Resulting from the prior knowledge of adult A. antennatus and the previous analysis, we released 20 000 particles in those areas on 1 and 15 July, 1 and 15 August, and 1 September.

The contribution of hydrodynamics to particle dispersions was explored through their horizontal and vertical displacements. The horizontal displacements drifted by particles (called drift distances) were the sum of the traveled distances between the daily recorded positions of particles. A Student two-sample test was computed to assess if the difference between those values according to a given factor (e.g., release zone, canyons, and hydrodynamics models) was significant.

In this study, the last position of the particles was attributed to the zone beneath the particle. The proportion of particles released from a release zone reaching a settlement zone was displayed in a connectivity matrix. Each cell represents the proportion of particles P_i,j from a zone i that has settled into a zone j:

$$\begin{array}{}\text{(4)}& P_{i,j}=N_{i,j}/N_i,\end{array}$$

where (i, j) is in the interval [1:10], N_i,j represents the number of particles settled in the zone j which has been released in zone i, and N_i is the number of particles released in zone i. Retention proportions are assumed to be the ratio of particles that remained in the zone where they were released (P_i,i with j=i) and appear on the diagonal of the matrix.

3 Results

Dispersal of particles in both configurations of ROMS had common general patterns and mostly diverged on the transport magnitude.

3.1 Lagrangian dispersal within the ROMS-Rutgers outputs

Releases in canyons generally made the drift distance larger and variable (Fig. 5). Lagrangian simulations carried out with the ROMS-Rutgers model transported particles over 27 km on average and up to 111 km for the longest trajectory. First, from the canyons, passive particles drifted 33 km, 6 km more than the average distance, while from open slopes, particles drifted 25 km, which was 3 km less than the average distance. Besides, the highest average drift distance was 36.3 km (±15.1 km) from the releases in the PC. Second, when released on the northern walls of the canyons, the particles drifted 29 km, which was 10 km less than the particles released on the southern walls. Last, the open slope 1 (OS 1) and OS 3 were the release zones with the shortest transports from which particles drifted 15.2 km (±5.5 km) and 15.7 km (±6.8 km), respectively.

Figure 5 Horizontal transport of passive particles (km) from the release zones in the ROMS-Rutgers (red bars) and ROMS-Agrif (white bars) simulations. Boxplots represent the minimum, 25 % quantile, mean, 75 % quantile, and maximum. Release zone codes are as indicated in Fig. 1.

Particle retentions in the release zones were high and particles from canyons generally seeded the nearby zones (see Fig. 6). The overall retention rate was averaged at 60 % (±33 %), because 42 % to 100 % of the particles stay in their release zones, except for the particles from the OS 2 (99 % of particles drifted into the PC). Particles released in the CCC can drift up to 28 % and 27 % in the southern zones of OS 2 and PC, respectively. The connectivity of CCC with PC is possible for particles having drifted approximately 50 km (±7 km) from 624 m (±92 m) depth. The zone OS 3 has received 56 % of particles from the PC and has retained 96 % of particles released inside. Last, we observe a drift direction opposite to the general southwestward direction of the main current in the south of our study area. Twenty-three percent of particles from the OS 4 drifted northward to the BC.

Figure 6 Dispersal rate (%) of passive particles between the zones of release (y axis) and settlement (x axis). The retention rates are on the diagonal. Release zone codes are as indicated in Fig. 1.

The vertical displacement of particles at the end of the 31 d simulation was relatively independent of the bottom depth at the release event (Fig. 7). Generally, more than 26 % of particles rise above the isobath of 500 m, around 30 % of particles drifted without having vertically displaced, and 5 % of particles go deeper than 800 m (see Fig. 7). More specifically, 12.8 % of particles released above the bottom between the isobaths of 700 and 800 m (hence, 700/800 m) arrived below 800 m. There are only 0.7 % and 4.8 % of bottom-released particles between the isobaths of 500 and 600 m and the isobaths of 600 and 700 m (hence, 500/600 and 600/700 m), respectively, which also reach such a depth. Moreover, the downward displacement is greater for particles released in canyons than on open slopes. From the three depth layers, 500/600, 600/700, and 700/800 m, as previously defined, 0.2 %, 1 %, and 11.1 % of particles on open slopes and 2 %, 13.5 %, and 19.7 % of particles on canyon slopes, respectively, go deeper than 800 m. Besides, the particles released in the middle layer (600/700 m) of the canyons were likely to be more dispersed vertically than from similar release conditions on the open slope (retention on the depth 27.3 % and 36.8 %, respectively).

Figure 7 Depth range (m) of particles after 31 d of drift. Particles were released between the following three depth ranges: 500–600 m (gray), 600–700 m (red), and 700–800 m (black) denoted by colored bars on the x axis.

In submarine canyons, the particles had a broader vertical displacement than over open slopes. Intensity and amplitudes in the vertical displacement of particles were dependent on the canyons and their walls, especially when they were released above both walls of the Blanes canyon (Fig. 8). On average, the particle vertical displacements in canyons or above the open slope are ascents of 28 and 24 m, respectively. Furthermore, the variability of the vertical displacement of particles in canyons (±154 m) is 2 times higher than on open slopes (±67.5 m). The drifts in PC have on average the highest ascent of 133 m. The PC is also the canyon with the widest vertical displacements where the maxima upward and downward displacements are 490 and 1163 m, respectively. From its southern wall, particles averagely rise 145 m, while on the northern wall there is an ascent of 58 m. Particles from the CCC averagely rose 44 m as did particles from the PC, but the vertical displacements were limited to a small range of 86 m in the CCC compared to 170 m in the PC. The drifts from BC have on average the furthest descent (127 m down). The drifts within BC went downward by 62 m when particles were above the northern wall and went upward by 36 m when particles drift above the southern wall. Contrary to the different days of release in the other canyons, the northern and southern walls of BC are where the temporal variability of the vertical displacements is the highest with a standard deviation of 53 and 34 m, respectively.

Figure 8 Particles vertical displacement. Vertical displacement according to particles drifting on the open slope or the canyon slope (a), drifting within the canyons (b), and released on the northern wall (North.) or southern wall (South.) of the canyons (c) in the ROMS-Rutgers and the ROMS-Agrif simulations. Release zone codes are as indicated in Fig. 1.

The influence of submarine canyons on the transport paths resulted in a common drift direction with small divergences in the trajectories (see Fig. 9). Indeed, particles southwardly followed the continental slope and spread differently according to the canyon in which they were released. To begin with, transport paths starting in the CCC and at one depth layer (500/600, 600/700, or 700/800 m) did not overlap with the ones starting at another depth layer, and they kept a linear direction when particles left the CCC. Second, the simulated drifts that began from the PC were affected by internal circulation within the canyon (i.e., recirculation of the particles between the walls), and transport paths intersected. Last, the transport paths initialized from the BC barely left the canyon area and overlapped with the trajectories set up in a different depth layer. The overlap or intersection of the transport paths may be a consequence from the vertical displacement of the particles (see Figs. 7 and 8) that were constrained by the horizontal and vertical circulation variability (see Figs. 5 and 8).

Figure 9 Estimated transport paths of passive particles during 31 d. Passive particles were released in the Cap de Creus canyon (a, d), Palamós canyon (b, e), Blanes canyon (c, f), and advected by the modeled hydrodynamics from ROMS-Rutgers (a–c) and ROMS-Agrif (d–f). Colors of the transport paths are associated with the depth ranges 500–600 m (gray), 600–700 m (red), and 700–800 m (black) on which particles were released.

3.2 Lagrangian dispersal within ROMS-Agrif outputs

The particles were transported further in the ROMS-Agrif configuration than in the ROMS-Rutgers configuration (Fig. 5), and particles from the canyon zones had traveled slightly further than particles from the open slopes. The general drift distances of particles are approximately 70.8 km (±26.3 km) and 44 km longer than drift in ROMS-Rutgers. Particles could drift up to the maximum distance of 197 km, detected for particles from the PC (see Fig. 5). Among the canyons, particles from CCC were transported the least distance: 47.1 km (±21 km). Oppositely, the released particles on PC averagely travel the furthest with drifts of 84.2 km (±24.9 km). Particles from canyon slope zones or open slope zones were transported 73.7 and 68.8 km, respectively, and even though the difference of drift distances is low (5 km), it is significantly different (p<0.05). Within the canyons, the particles from the northern walls of the CCC and PC drift 66 km, which was 13 km less than the particles from the southern walls. In the BC, the tendency was for particles released on the northern wall to travel 84 km, or 10 km more than particles from the southern wall.

Southern transport of particles was on a broader scope than in ROMS-Rutgers (see Fig. 6). After a simulation period of 31 d, the retention rate by the release zones is 1 % to 43 % of particles (Fig. 6), making an average of 18 % (±14 %) of particles. Like in ROMS-Rutgers, the open slope OS 2 retains the least of particles (1 %). Consequently, a small proportion of 29 % of particles ends up in the PC, while the majority (64 %) connects with open slope OS 3. The particles from the CCC also have drifted in more zones than any other release places. They were transported in four different zones, which were OS 1 at 6 %, OS 2 at 9 %, PC at 45 %, and OS 3 at 19 %. Particularly, 45 % of particles connect to the PC and those particles are characterized by a drift of 70.5 km (±9 km) from 575 m (±77 m). A lower rate of dispersal is detected, as well, from PC to BC. On average, 25 % of particles drift approximately 86.6 km (±14.9 km) from 600 m (±82 m) to reach BC.

The pattern of vertical particle displacements in ROMS-Agrif was similar among the different release depths and topographic structures of the NW Mediterranean Sea (see Fig. 7). Overall, 37.3 % of particles reaches the layer above 500 m depth, while only 0.8 % of particles go deeper than 800 m. Particles released over the upper bottom depths (500/600 and 600/700 m) presented comparable but still significantly different (p<0.05) vertical displacements even if they were released in canyons or on open slopes. Around 41 % and 43 % of particles drifted within the same depth layer as their releases and 44 % to 48 % ascended in upper depths. On the contrary and overall, from releases between isobaths 700/800 m, a majority of the particles (65 %) went upward, a high proportion of particles (6.1 %) went down below 800 m, and fewer particles were transported in the depth layer of their release (retention in the deeper layer of 29 %).

The average vertical displacement of particles was a small ascent in the water column (see Fig. 8), with canyons being the places of higher ascents. On average, particles rose 15 m (±53 m) when released on open slopes, and twice that, 34 m (±65 m), when released in canyons. Among the different canyons, particles released in the PC were vertically displaced the least (23 m on average, or half the rise of particles from the CCC and BC), but the widest vertical displacements of 352 m downwards and 519 m upwards were from PC. We furthermore notice that the average ascent and descent of particles were the highest in BC (74 m) and in PC (44 m), respectively, which are reversed directions of movement compared to the results in ROMS-Rutgers. Lastly, particles drifting above the northern wall of all the canyons significantly (p<0.05) rise higher than for particles above the southern wall. In CCC, the rise of particles was 2.7 times higher on the northern wall (71 m) than on the southern wall, which represents the highest change in displacement between the walls. In contrast, in PC, the rise of particles reaches 68 m on the northern wall, which was 4 m less than for particles on the other wall. Over the different release dates, the southern wall of PC is where the vertical movement of particles temporally varies the most (standard deviation of 39 m).

With the simulation of passive drifts in ROMS-Agrif, patterns of particle trajectory were similar to the ones in ROMS-Rutgers, especially the southward drift direction and the spread of particles from the submarine canyons (see Fig. 9). Nonetheless, the length of the transport paths, as demonstrated by the drift distance values (see Fig. 5), was longer than in ROMS-Rutgers, and the spread of trajectories was wider. The spread was particularly amplified when particles drifted in the Palamós and Blanes canyons. For example, transport paths from the PC disclosed the advection of particles taken in two circular flows near this canyon, and the span of the transport paths from the BC widened once the particles left the canyon. Like for the simulations in ROMS-Rutgers, the variability in the transport path originated in the depth position of the particles (see Figs. 7 and 8). However, the spread of the trajectory was intensified as a consequence of the current intensity estimated in the ROMS-Agrif (see Fig. 4).

4 Discussions

Two hydrodynamic models with different implementations allowed for testing the sensitivity of particle dispersals to the Mediterranean Sea thermohaline circulation. The dispersal rates were low and similar to observations or simulations of deep larvae from other marine ecosystems like in the Baltic Sea and in another basin of the Mediterranean Sea (Corell et al., 2012; Palmas et al., 2017).

4.1 Deep-sea connectivity between submarine canyons

For the first time in the NW Mediterranean Sea, the connectivity between submarine canyons has been demonstrated through the Lagrangian transport of deep particles. The southwestward circulation having important amplitude in upper layers also predominated in the bottom drift simulation and explained the particle trajectories from the northern canyons to the southern canyons. On a local scale (i.e., canyon area), submarine canyon areas received drifters from upstream release and seeded the nearest downstream open slopes. In fact, this result is related to the flows of the NW Mediterranean Sea following the right-hand side canyon walls (Flexas et al., 2008; Palanques et al., 2005; Granata et al., 1999; Jorda et al., 2013). With this, our drift simulations supported one of the submarine canyon functions, which consisted of seeding particles in neighboring zones.

Nonetheless, each submarine canyon has a singular topography and a different angle of exposure to the main circulation, resulting in different patterns of particle dispersals. In the horizontal dimension and inside each submarine canyon, particle drifts were different because they did not begin from the same depth, the same day, or with the same exposure to the current (i.e., releases on upstream or downstream canyon walls) even though the general canyon influence led to common north–south drift patterns. Those differences also appeared in analyses of sediment flux in the Mediterranean Sea (Durrieu de Madron et al., 1999; Lopez-Fernandez et al., 2013) and in simulations of passive particle drifts in different water layers (Ross et al., 2016). In the vertical dimension, our results reflected the uncorrelated currents occurring inside a canyon (Jorda et al., 2013) that could be one of the origins of the drift differences between canyons. For example, the amplitude of vertical drifts in Palamós canyon was greater than vertical drifts in the other canyons. The uncorrelated currents could also explain the general upward displacement of particles that were released on the upstream canyon walls. On the contrary, in shallow waters, previous flow studies on those canyon sides estimated a downward flux (Ahumada-Sempoal et al., 2015; Flexas et al., 2008; Granata et al., 1999).

Even though the main trajectory patterns persisted, canyon influence on dispersals changed over the summer releases. In Kool et al. (2015), the main deep circulation pattern over 3 years was consistent enough to not strongly impact the dispersals. This finding also applied for our results when the range of releases was over 2 months. Temporal variability around the main patterns was likely caused by the dynamic of the internal circulation over time inside the canyons (Ahumada-Sempoal et al., 2013; Palanques et al., 2005; Jorda et al., 2013). Nonetheless, further simulations of drifts need to be done to identify the minimum frequency of release to observe significant changes in the drift variability.

4.2 Sensitivity of the passive transport to deep-sea circulations

In our study, the large amplitude of dispersals in ROMS-Agrif was the main origin of differences between the drifts. Particle retention was 3.3 times higher with the simulation carried out with ROMS-Rutgers than with ROMS-Agrif. Whether the particles were released on canyon slopes or on open slopes did not change the fact that distance drifts were 2 times longer in simulations from ROMS-Agrif. The largest difference in amplitude of drifts (an average difference of 64.2 km) was simulated from the open slope between Palamós and Blanes canyons. There, in the ROMS-Rutgers bathymetric grid, the presence of a submarine valley in the slope (∼41.5^∘ N, ∼3.25^∘ E) lowered the dispersal rates, and acted as a wall that generated distinct dispersal lengths depending on the side on which particles were released. In contrast, on the open slope between the Cap de Creus and Palamós canyons, the continental slope did not have specific topographic structure, and in both models the zone retention was low.

Thus, small topographic differences (e.g., the submarine valley) indirectly influenced the particle transport by conditioning the flow circulation when the hydrodynamic models were run. The main cause of topographic differences was related to the resolution of the bathymetric data sets chosen by the modelers. Indeed, an adequate bathymetric data set would approximate better the mesoscale structures of the water circulation (Gula et al., 2015). Then, it made the simulated drifts from the ROMS-Rutgers more appropriate for the deep passive drifts because the bathymetric data set was more precise (i.e., horizontal resolution of 100 m instead of 3 km) even though it has a lower horizontal resolution than the ROMS-Agrif. Another major difference within the drifts in the two hydrodynamic models was the northward drifts from Blanes canyon and from the open slope at its southern end (OS4) in ROMS-Rutgers. The northward particle drift is a consequence of an anticyclonic flow generated at the mouth of Blanes canyon that also affects the local deep water circulation (Ahumada-Sempoal et al., 2013; Jorda et al., 2013). The southwestward flow in ROMS-Agrif may be the consequence of the bathymetric data set used before the run of the hydrodynamic model and the hydrodynamic model parameterization of the water profile discretization. In ROMS models, the water is discretized by terrain-following sigma layers that generally have high resolution near the surface and relax the circulation variability near the bottom.

The generated drift uncertainties gave arguments for choosing one hydrodynamic model over another, above all when drifts are in deep water. In the absence of more hydrodynamic models, we encourage the deep-sea Lagrangian modeler to cross-validate the Lagrangian drifts with the existing literature. Our work evidences the importance of two major sources of uncertainty in the particle drifts from the hydrodynamic models: bathymetry and model discretization (i.e., choice of S-coordinate parameters). These two factors are shown to be particularly important in the zones near the canyons (all the water column) and open-slope areas. This is an important outcome for modeling studies. In this sense, this study shows the areas where uncertainty grows and so helps to evaluate when and where the choice of the hydrodynamic models will have a large impact on the drift results.

4.3 Relation between deep-sea circulation and larval cycle of A. antennatus

The association of eggs and larvae of Aristeus antennatus with passive particles comprised the first approach to their drift. Our study revealed interesting dispersal features which could be related to the ecology of A. antennatus. First, few particles in canyons were transported near 1000 m depth, where peaks of juvenile abundances have been reported at the end of fall (Sardà et al., 1997; D'Onghia et al., 2009; Sardà and Cartes, 1997). It would imply that the late larvae are helped or constrained by the vertical circulation to settle in the deep zones. Second, the low dispersal rates from the canyons highlighted that inactive larvae may be retained by those topographic structures. In fact, the morphology of the canyons facilitates the retention of pelagic particles for a few days (Ahumada-Sempoal et al., 2015; Rojas et al., 2014). In that case, it would mean that subpopulations of A. antennatus strongly depend on their own. Then, the surface of management plans should be structured to local shrimp fishing grounds, like the one in Palamós canyon (Boletín Oficial del Estado, 2013, 2018). Furthermore, because the submarine valley is localized between two areas of high fishing effort (Palamós and Blanes), there is an interest in its potential role for marine populations. Literature categorizes submarine valley effects as canyon influences, but, due to the exposure of the submarine valley to the current, an evaluation of whether the valley deviates or recirculates the flow is needed.

In this study, the duration of simulation corresponded to the first PLD approximation of Aristeus antennatus larvae generalized from penaeid larvae. Nonetheless, further studies need to explore the variability of the estimated dispersal features with longer PLDs. In the framework of the A. antennatus case study, the PLD length can be expanded in accordance with a plausible time range of 6 months, i.e., time between the beginning of the spawning peak in July and the period of sampled early juveniles (Sardà et al., 1997; D'Onghia et al., 2009; Sardà and Cartes, 1997). The effect of PLD length can also be estimated with water temperature changes during the larval drift, as they are likely to affect the PLD when A. antennatus larvae reach the warm surface waters (O'Connor et al., 2007). The larval behavior needs to be taken into account, because it can significantly affect the vertical position of decapod larvae in the water column and may influence the larval dispersal (Cowen and Sponaugle, 2009; Levin, 2006; Queiroga and Blanton, 2004). Few studies on A. antennatus larvae sampled in open ocean have revealed their presence in the surface water layer (Carbonell et al., 2010; Carretón et al., 2019; Heldt, 1955; Seridji, 1971; Torres et al., 2013). Sardà et al. (2004) and Palmas et al. (2017) assumed that positive buoyancy of eggs partly underlined the presence of individuals in the shallowest water layer. The buoyancy mechanism needs to be analyzed through sensitivity tests in order to compensate for the lack of accurate knowledge about it (Hilário et al., 2015; Ross et al., 2016).

5 Conclusions

Passive drifts of particles simulated by a Lagrangian transport model coupled with two hydrodynamic models (ROMS-Rutgers and ROMS-Agrif) were compared. The particle releases used characteristics of deep-sea red shrimp, A. antennatus, spawners and scarce knowledge of their free-living period (i.e., eggs and larvae). Then, drift simulations began in a highly disturbed circulation induced by the submarine canyons of the NW Mediterranean Sea. In general, horizontal dispersion of deep-sea passive particles was relatively short during the 31 d simulation, but their transports could establish links between the submarine canyons. Like numerous other studies, our work concludes that marine topographic features (i.e., the deep submarine canyons and submarine valleys) favored particle reception and retention. It also supports that the amplitude and variability of vertical displacement (upward and downward in the water column) are higher inside submarine canyons than on other parts of the continental slope. Besides, the variability in the particle vertical displacement was dependent on the three-dimensional distribution of particle releases in a submarine canyon. Variability in particle dispersions was also a consequence of differences between canyon topography and their exposure to the currents. This sensitivity to the topography highlighted the importance of the modeler choice before running hydrodynamic models. Having a finer bathymetry data input and a finer vertical resolution near the bottom (i.e., like in ROMS-Rutgers) is important for estimating the circulation and the passive dispersions with more accuracy, even though the hydrodynamic model has a coarse horizontal resolution (i.e., 2 km vs. 1.2 km). Therefore, the ROMS-Rutgers model was approved by better modeling the bottom drifts, and it will be used for future studies on the benthic deep-sea species A. antennatus. Numerical simulations of the larval drift are one of the best approaches that allowed for approximating the unknown dispersal paths of A. antennatus in the NW Mediterranean Sea. This study is the first step of a work conducted with the final focus to improve fishery management of a deep-sea species. However, observation of the pelagic larval traits, validation of the dispersals, and simulation including biological behaviors (vertical migration, buoyancy, and swimming abilities) are required before involving the results in the management strategies.