# Sensitivity and robustness of larval connectivity diagnostics obtained from Lagrangian Flow Networks

Pedro Monroy, Vincent Rossi, Enrico Ser-Giacomi, Cristóbal López, Emilio Hernández-García

► To cite this version:

Pedro Monroy, Vincent Rossi, Enrico Ser-Giacomi, Cristóbal López, Emilio Hernández-García. Sensitivity and robustness of larval connectivity diagnostics obtained from Lagrangian Flow Networks. ICES Journal of Marine Science, 2017, 10.1093/icesjms/fsw235. hal-01731177

HAL Id: hal-01731177

https://amu.hal.science/hal-01731177v1

Submitted on 24 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sensitivity and robustness of larval connectivity diagnostics obtained from Lagrangian Flow Networks

Pedro Monroy

IFISC (CSIC-UIB), Palma de Mallorca,
07122, Spain

pmonroy@ifisc.uib-csic.es; tel: +34971259881; fax: +34971173248

Vincent Rossi

IFISC (CSIC-UIB), Palma de Mallorca,
07122, Spain

vincent@ifisc.uib-csic.es

Enrico Ser-Giacomi

Ecole Normale Supérieure, PSL Research University, CNRS, Inserm,
Institut de Biologie del École Normale Supérieure (IBENS), F-75005
Paris, France

enrico.sergiacomi@gmail.com

Cristóbal López

IFISC (CSIC-UIB), Palma de Mallorca,
07122, Spain

clopez@ifisc.uib-csic.es

Emilio Hernández-García

IFISC (CSIC-UIB), Palma de Mallorca,
07122, Spain

emilio@ifisc.uib-csic.es
**Abstract**

Lagrangian Flow Network (LFN) is a modeling framework in which geographical sub-areas of the ocean are represented as nodes in a network and are interconnected by links representing the transport of propagules (eggs and larvae) by currents. Here we assess the sensitivity and robustness of four connectivity metrics derived from LFN that measure retention and exchange processes, thus providing a systematic characterization of propagule dispersal. The most relevant parameters are tested over large ranges: the density of released particles, the node size (spatial-scales of discretization), the Pelagic Larval Duration (PLD) and the modality of spawning. We find a threshold for the number of particles per node that guarantees reliable values for most of the metrics examined, independently of node size. For the region and dates considered in this study this threshold is about 100 particles per node. We also find that the size of network nodes has a non-trivial influence on the spatial variability of both exchange and retention metrics. Although the spatio-temporal fluctuations of the circulation affects larval transport in a complex and unpredictable manner, our analyses evidence how specific biological parametrization impact the robustness of connectivity diagnostics. Connectivity estimates for long-PLD organisms are more robust against biological uncertainties than for short-PLD ones. Furthermore, our model suggests that for mass-spawners that release propagules over short periods (≈ 2 to 10 days), daily release must be used to properly consider connectivity fluctuations. In contrast, average connectivity estimates for species that spawn repeatedly over longer duration (a few weeks to a few months) remain robust even when using longer periodicity (5 to 10 days). Our results have implications to design connectivity experiments with particle-tracking models and to evaluate the reliability of their results.

**Keywords:** Bio-physical model; particle-tracking model; Connectivity; Lagrangian Flow Network; Larval dispersion and transport; Sensitivity analysis; Robustness
1

# Introduction

2

3

4

5

6

7

8

9

10

11

12

Dispersal, the tendency for organisms to move, settle and ultimately reproduce away from their native habitats and parents, influences many evolutionary and ecological processes, from speciation and species persistence (long-term survival) to biological invasions and the resilience of ecosystems to external perturbations (Duputie & Massol 2013). Dispersal and the tightly related concept of “connectivity” have been shown to favour the access to feeding and nursery grounds (Mbaye et al. 2015), to influence demographic rates of sub-populations (Palumbi 2003), to control the coexistence of genotypes within species (Lowe & Allendorf 2010), to affect population structures and dynamics (Cowen & Sponaugle (2009); M. Hidalgo, pers. comm.), to regulate biodiversity (Clayton et al. 2013; Jones et al. 2007) and to plan conservation strategies (e.g. Guizien et al. 2012; Garavelli et al. 2014; Rossi et al. 2014; Dubois et al. 2016).

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

Studying connectivity in marine ecosystems, as in terrestrial ones (Carr et al. 2003), is thus crucial. Yet, the highly variable properties of seawater and its dynamical character have profound effects on the connectivity of marine systems. Indeed, because the life cycle of most marine organisms encompass an embryo phase (eggs and larvae), and considering that juveniles and adults tend to be rather territorial, one of the most important process of dispersion is the transport of propagules by ocean currents over large distances (Cowen & Sponaugle 2009). The marine realm is an heterogeneous medium in perpetual movement with non-linear processes occurring at various scales, from diffusive dynamics, eddy-driven mixing, up to basin-scale advective currents (Rossi et al. 2013; Clayton et al. 2013). Over the last decades, oceanographers recognized the prominence of mesoscale and sub-mesoscale structures (whose scales span approximatively over 1-100 km and 1 day to several months) in oceanic flows (Rossi et al. 2008; Capet et al. 2008; Mancho et al. 2008; Mahadevan 2014). These studies revealed the presence of barriers (which prevent the dispersion of propagules) and of preferential corridors of transport (favouring efficient connections), having important consequences for the dispersal and retention of non-swimming larval stages (Harrison et al. 2013; Vaz et al. 2013; Rossi et al. 2014; Dubois et al. 2016). Even for life-stages that are actively mobile, the oceanic circulation governs the spatial distribution of environmental variables and planktonic food (d'Ortenzio & d'Alcalá 2009), thus playing a crucial role in controlling the location, timing and success of spawning, dispersal and settlements (Mariani et al. 2010; Reglero et al. 2012).

33

34

35

36

37

38

Even though marine connectivity has been studied from several approaches, it remains difficult and costly to evaluate (Calò et al. 2013). Sparse observations of larval distribution do not provide enough information to infer broad-scale connectivity. Genetic methods provide estimate of connectivity, usually at long -evolutionary- time-scales, but costly sampling tends to allow only for poor spatial coverage, preventing the full description of spatial connectivity patterns. In this respect, biophysical models are relatively affordable
and extremely powerful to study and forecast dispersal patterns at different temporal and spatial scales. As such, a significant part of connectivity research employs numerical bio-physical models, especially particle-tracking (or Lagrangian) models, to characterize transport and dispersion in the ocean (e.g. Huggett et al. 2003; Mitarai et al. 2009; Treml et al. 2012; Andrello et al. 2013; Wood et al. 2014; Kool & Nichol 2015). Despite the widespread use of these modelling approaches, there is still a poor knowledge of the key biological traits (i.e. spawning areas, larval duration, larval behavior and mortality rates) from which most model parameters are derived (Calò et al. 2013). In addition, it is not clear how different parameter choices could affect connectivity estimates and their variability. Last, in spite of the recognized influence that ocean currents have on the connectivity of marine ecosystems, these processes remain difficult to evaluate due to their inherent chaotic, unpredictable and multi-scale properties. For instance, Huggett et al. (2003) reported that the transport patterns of anchovy eggs are most sensitive to changes both in areas and in months of release (along with the corresponding changes in circulation) than to variations in biological parameters. Other species-specific studies found that the temporal variability in circulation is crucial in modulating dispersal (Qian et al. 2015), while others found that biological processes crucially affect modelled connectivity estimates (e.g. Parada et al. 2003; Paris et al. 2007). These examples evidence the complexity of simulating tightly coupled bio-physical processes and the importance of parametrization. It calls for further quantitative and systematic characterization of the intrinsic functioning and variability of these models.

With the growing use of connectivity models to estimate connectivity of a given organism, scientists realized the importance of assessing their robustness and sensitivity. Most existing studies focused on the evaluation of uncertainties associated with the underlying hydrodynamical model (e.g. Putman & He 2013) or of those related to changes in any bio-physical parameters (such as larval mortality, vertical migration, depth of dispersion, density and frequency of particle releases) but considering a specific species in a given oceanic region (e.g. Huggett et al. 2003; Paris et al. 2007; Queiroga et al. 2007; Huret et al. 2010; Treml et al. 2012; Cuif et al. 2014; Garavelli et al. 2014; Kough & Paris 2015). This implicitly means that those results would only be valid for the species and region of interest. In contrast, our objectives are to explore the sensitivity of connectivity proxies to the commonly used parameters over large ranges of values without restricting our analyses to a target species. We also focus on an extended and fully covered region (Mediterranean basin) examining systematically the influence of different discretization-scales on the spatio-temporal sensitivity of the aforementioned diagnostics.

To our knowledge, Simons et al. (2013) is the only study analysing systematically the sensitivity of larval transport predictions to large ranges of three universal input parameters: the number of particles released, the particle release depth and the particle tracking time. They showed that larval transport predictions are sensitive to changes in
all three input parameters and are dependent on the local mixing strength. While Simons
et al. (2013) is the first sensitivity analysis of larval transport of this broad type, it solely
focused on a dozen of near-shore sites over a relatively small oceanic region (the Southern
California Bight). Moreover, this investigation is restricted to two-dimensional vertically-
integrated Particle Density Distributions (PDDs) and to the 2-dimensional grid-based
descriptors of connectivity directly related to them (such as raw connectivity matrices,
dispersal kernels, probability density functions and particle concentrations).

We aim here at studying the general robustness and sensitivity of four connectivity proxies derived from Lagrangian Flow Networks (LFNs) (Dubois et al. 2016). The LFNs modelling framework (Ser-Giacomi et al. 2015a) has demonstrated great effectiveness and flexibility in characterizing transport phenomena in geophysical flows due to its solid theoretical background relating network measures with dynamical properties of the flow (Froyland et al. 2007, 2012; Ser-Giacomi et al. 2015a). This paradigm has been applied to identify hydrodynamical provinces which organize larval dispersion (Rossi et al. 2014) and to interpret local connectivity measures in the context of population structuring (Dubois et al. 2016) allowing in such way a multi-scale characterization of the dispersion of planktonic organisms. The LFN framework permits also unveiling the main avenues of transport through path analysis (Ser-Giacomi et al. 2015b; Jönsson & Watson 2016). While this approach has been already applied in different contexts, the global robustness and sensitivity of metrics derived from LFNs have not been quantitatively assessed yet.

The four complementary connectivity proxies tested here measure the retention and exchange of larvae and are calculated by manipulating elements of the classical connectivity matrix and its binary equivalent (see also sect. *Materials and Methods* and Dubois *et al.* 2016). Through their integrated analysis, they provide more information than, and are not directly comparable to, classical connectivity diagnostics mentioned previously (Dubois *et al.* 2016). In addition, we extend the results of Simons *et al.* (2013) by focusing on all sites of a larger and more heterogeneous oceanic region, namely the Mediterranean basin, which cover a large range of flow patterns and mixing intensities.

In this study we assess the sensitivity and robustness of connectivity proxies derived
from LFNs over the Mediterranean basin. After presenting briefly the LFNs methodol-
ogy and reporting the statistical methods used (sect. *Materials and Methods*), we explore
the sensitivity of the four metrics to numerical parameters, namely the density of initial
particles and the node-size (sect. *Numerical parameters*). Then, we evaluate the impact
of biological parameters such as the Pelagic Larval Duration (PLD), simulated by the
period of tracking, and the spawning dates, simulated by the starting times (sect. *Bi-
ological parameters*). Indeed, from the perspective of marine ecologists, there are high
uncertainties in the appraisal, based on field observations or laboratory experiments, of
the key life history traits controlling connectivity. These include estimates of the PLD
(MacPherson & Raventós 2006), of the modality of spawning (Kough & Paris 2015) and
of the depth of dispersion (Corell et al. 2012; Cuif et al. 2014). Our results are interpreted and discussed against existing literature (sect. *Discussion and conclusions*) and we finally conclude with recommendations for future connectivity studies using particle tracking models and with foreseen impacts of such parameter choices on connectivity metrics.

## Materials and Methods

### The Lagrangian Flow Network framework

We study larval dispersion with the methodology of Lagrangian Flow Networks (LFN) (Rossi et al. 2014; Ser-Giacomi et al. 2015a). It represents fluid transport among different oceanic regions as a network of nodes interconnected by directed and weighted links. Each node corresponds to a geographical sub-area of the oceanic surface and can be seen as a discrete habitat patch. Each link symbolizes an effective mass transport driven by ocean currents between any two of these sub-areas during a given time interval, representing larval transport pathways. This framework is ideal to study the structure of marine populations which are commonly organized as discrete subpopulations whose respective persistence depends on self-replenishment and asymmetric exchanges of individuals (here assumed to be mainly mediated by larval transport) among fragmented habitats (nodes) (Kool et al. 2013; Dubois et al. 2016).

Applying the LFNs to a given oceanic area requires the computation of Lagrangian trajectories of fluid particles (simulating passively drifting larvae) driven by an eddy-resolving velocity field of the region. Two-dimensional trajectories are obtained at a fixed depth through the integration of the velocity field, linearly interpolated in time and space, using a fourth order Runge-Kutta algorithm (Ser-Giacomi et al. 2015a) with integration time step 1 day. All outputs presented in this paper used the daily horizontal flow field generated by Copernicus operational model implemented in the entire Mediterranean Sea at 1/16° horizontal resolution over 26 years (1987-2012) (Oddo et al. 2009). We retained the surface layer at 12 m because it is arguably the most turbulent flow field (ocean interior circulation becomes more sluggish) that is not directly influenced by wind drag (as is the very first surface layer). It means that our results are on the conservative side: computation at deeper layers would return similar or even more robust metrics than those computed at 12 m because of the reduced turbulence and variability. Integrating particles in three-dimensions and the very few trajectories that end up on land before the integration ends (“beaching”) have barely detectable impacts on our results, as discussed in Rossi et al. (2014), Ser-Giacomi et al. (2015a) and in section *Discussion and conclusions*.
Besides the oceanic domain of interest, the selected depth for deploying the particles and the corresponding velocity field, the main parameters required by the LFN methodology are the node size (i.e. scale of spatial discretization), the area density of particles released within every node, the spawning time(s) at which particles are released and the duration of tracking which simulates the Pelagic Larval Duration (PLD, i.e. time larvae drift with ocean currents). This paper analyses the sensitivity of four connectivity estimates to those parameters varying over large ranges of variability, within the limits imposed by our settings (table 1).

The construction of LFNs consists in four steps:

1. Divide the region in oceanic boxes which are the nodes of the network (the size of these boxes is in general larger than the grid used by the hydrodynamic model to compute velocities) and then initialize evenly the particles inside each node.

2. Compute the Lagrangian trajectories for all particles in each node and record the node of their final positions.

3. Calculate the connectivity matrix (or adjacency matrix) **P** of the LFN and its binary equivalent **L**. The element $P_{i,j}$ of the connectivity matrix is defined as the probability that a particle released in node *i* settles at node *j* after the tracking time (PLD). It can be estimated by the number of particles exchanged from node *i* to node *j* divided by initial number of particles in node *i*. The elements of the binary connectivity matrix are $L_{i,j} = 0$ if $P_{i,j} = 0$ (i.e. no larval flow) and $L_{i,j} = 1$ when $P_{i,j} \neq 0$ (Ser-Giacomi et al. 2015a).

4. Perform various and complementary post-processing of both matrices (**P** and **L**) to obtain relevant metrics and study connectivity.

The framework is tunable to any target species through the use of specific parameters that match the traits of this organism. For instance, the dispersal of hake's larvae are being studied with the LFN by simulating successive spawning events during autumn using a PLD of 40 days and the currents at ~ 100 m (M. Hidalgo, pers. comm.). Here we used LFN in a generic configuration for an ecosystem approach to connectivity using ecologically-relevant values of parameters (Dubois et al. 2016). Building upon the bibliography review of Dubois et al. (2016), we explore large range of parameters to derive generic relationships and to relate them to the physical properties of the underlying flow field. We also assume that local larval release, larval mortality and success of recruitment are homogeneous in space in order to test the robustness of connectivity estimates for the most generic biological system.<fnref n="1" /> This original perspective provides insights about the

<footnote n="1">¹Note that, within the LFN framework, these neglected biological processes can be easily taken into account by manipulating “a-posteriori” the original connectivity matrix. For instance, one can incorporate heterogeneous larval production, larval mortality and success of recruitment by multiplying the generic matrix by diagonal matrices containing the node-specific magnitude of the process to model.</footnote>
intrinsic sensitivity and help choosing adequately parameters in our settings as well as in
any other off-line Lagrangian models.

This study focuses on local (node-scale) features: the fourth step before-mentioned
consists in the computation of four complementary metrics inspired from concepts of
Population Dynamics. Each connectivity matrix is first manipulated to examine, for
each node, the number and fluxes of all links emanating from it and entering it. Then,
through a combination of these basic variables, four indices of retention and exchange of
larvae are computed using the formulae reported in Dubois et al. (2016):

• Local Retention ($LR_i$, computed from **P**) is the ratio of locally produced settlement to local larval release, i.e. the proportion of locally retained particles in each node<fnref n="2" />.

• Self-Recruitment ($SR$, computed from **P**) is the ratio of locally produced settlement to the overall settlement, i.e. the ratio of retained particles in each box to the total incoming particles from all origins (including those produced locally).

• Source-Sink degree ($SS_d$, computed from **L**) evaluate the source/sink character in terms of the diversity of origin and destination sites.

• Source-Sink strength ($SS_s$, computed from **P**) evaluate the source/sink character in terms of the total fluxes of larvae (relative importance of larval export versus import).

These node-scale metrics are all restricted to values between 0 and 1. Elevated retention is represented by *LR* and *SR* approaching 1. The negative (positive) deviations of *SS*<sub>d</sub> and *SS*<sub>s</sub> from 0.5 allow the quick appraisal of sources (sinks).

Another relevant quantity attributed to each node is the surface proportion covered
by land (land ratio). It distinguishes “fully oceanic” nodes from coastal ones. As an
example, the land ratio of all nodes at 12 m depth are shown in Fig. 1.

Sensitivity tests and statistical methods

Sensitivity tests are performed for all tunable parameters which are distinguished in two
types. The numerical parameters refer to the node size and the density of particles
initially released in each node. We want to understand how those intrinsic parameters
and their variations affect the computation of connectivity metrics with the LFNs. With
respect to the biological parameters, we will focus on perhaps the most universal life
history traits used in dispersal model: the PLD, the date and frequency of spawning.
Note that our analyses are restricted to the parameters commonly employed in off-line
particle tracking models, while we did not explore the resolution and parametrization

<footnote n="2">²In the following we neglect the sub-index *i* which refers to any arbitrary node.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Grid used to construct the LFN of the Mediterranean Sea with node size of 1/4° and labelling the nodes according to its land ratio at 12m. The discrete colour-scale represents the value of the land ratio and the black line is the 12m isobath from the ETOPO1 dataset. The interior of the red rectangle is the region used in all analyses except the one in Fig. 2</figcaption></figure>

of the velocity field and corresponding ocean model (see also sect. *Other factors not considered*).

Quantitative comparisons among pairs of different experiments (e.g. 2-dimensional fields of any of our four connectivity metrics generated with two different set of parameters) are made by computing Pearson correlation coefficient *r*. It provides an absolute scale to measure similarity of spatial distributions of connectivity metrics obtained varying one parameter: In general $0 \le |r| \le 1$, with $|r| = 1$ indicating a perfect linear correlation, and $r = 0$ absence of linear correlation. Since *r* compares the value of two discrete fields node by node, it is highly sensitive to changes in spatial patterns. Note that this can be a disadvantage when one wants to focus on temporal variability with small nodes because when pattern suffers a small translation or rotation in space while their values remain practically the same across two experiments, the resulting *r* will be small while the patterns are visually identical.

It is worth noting that Pearson correlation measures linear relationship, so that its use may be questionable if the variables being compared follow a different relationship. It is also sensitive to outliers that might lead to large effects on *r* value. In order to check linearity and presence of outliers, we visually inspect the cloud of points on the scatter-plots of corresponding metric values in each node (e.g. Fig. 2). Additionally we
<table><caption>Table 1: Summary of experiments performed to assess the sensitivity and robustness to parameter changes of the four connectivity metrics derived from LFNs. Trajectories were computed using a time step of one day and the native resolution of the hydrodynamical model 1/16°.</caption><thead><tr><th>Sensitivity test</th><th>Parameter</th><th>Range value</th></tr></thead><tbody><tr><td rowspan="5">Node size and particle density (total: 1729 experiments)</td><td>Starting dates</td><td>From 1 September 2010 to 30 November 2010, every 5 days</td></tr><tr><td>PLD</td><td>60 days</td></tr><tr><td>Particle densities (km<sup>-2</sup>)</td><td>0.026, 0.052, 0.259, 0.518, 1.037, 1.556, 2.07, 2.588, 3.109, 3.630, 4.159, 4.674 and 5.184</td></tr><tr><td>Node sizes</td><td>1°, 3/4°, 1/2°, 1/4°, 1/8°, 1/10° and 1/12°</td></tr><tr><td>Depth</td><td>-12m</td></tr><tr><td rowspan="5">Starting time and spawning frequency (total: 364 experiments)</td><td>Starting dates</td><td>From 1 September 2010 to 30 November 2010, every day</td></tr><tr><td>PLDs</td><td>30 and 60 days</td></tr><tr><td>Particle density (km<sup>-2</sup>)</td><td>0.518 (i.e. 1600 and 100 particles per node of 1/2° and 1/8°, respectively)</td></tr><tr><td>Node sizes</td><td>1/2° and 1/8°</td></tr><tr><td>Depth</td><td>-12m</td></tr><tr><td rowspan="5">Uncertainties in PLD estimates (total: 1200 experiments)</td><td>Starting dates</td><td>1 September and 1 October each year over 2001-2010</td></tr><tr><td>PLDs</td><td>from 5 to 65 days, with daily increments</td></tr><tr><td>Particle density (km<sup>-2</sup>)</td><td>0.518 (i.e. 1600 particles per node of 1/2°)</td></tr><tr><td>Node size</td><td>1/2°</td></tr><tr><td>Depth</td><td>-12m</td></tr></tbody></table>

tested the non-parametric Spearman correlation (not shown) which is known to be more robust against outliers. We verified that our results are basically equivalent to the ones obtained with the Pearson correlation.

Over the range of parameter values upon which sensitivity is tested, a reference simulation should be chosen considering the optimal values that would return the “best” connectivity estimates. Then sensitivity experiments are generated by changing one parameter at a time and are finally compared against the reference. For instance, to explore the robustness of connectivity metrics to the density of particles, we chose the experiment using the largest number of particles as the reference simulation, as in Simons et al. (2013). Then, simulations obtained with decreasing number of particles are compared with the reference one.

Overall, the ensemble of numerical experiments performed are summarized in table 1 and are described in the following sections.
## Numerical parameters

We explore the sensitivity of connectivity diagnostics with respect to the density of particles initially released to determine the best compromise density that ensure robust estimates and reasonable computational cost. Simulations were conducted releasing particles every 5 days from 1 September 2010 to 30 November 2010 (a total of 19 starting dates) with a set of 13 different particle densities (spanning 0.02 up to 5.17 $km^{-2}$) and 7 different node-sizes (ranging $1/12^\circ$ - 1°). Details are reported in table 1. We use a PLD of 60 days for every simulation because it is considered as the longest integration time for which the 2-dimensional assumption (neglecting vertical movement) remains accurate (Rossi et al. 2014). Since longer PLD tend to favour long-distance dispersal (Shanks 2009), especially for passive larvae, it can also be considered as the “worst-case” scenario (i.e. results for shorter PLDs will return higher correlation coefficients, i.e. will be less affected by particle density).

We compare each experiment to the reference configuration with the largest density of particles ($5.17 km^{-2}$) computing Pearson correlation for the four connectivity metrics. For every pairwise correlation corresponding to each specific pair of node-size and particle density, a temporal average was performed among the 19 different starting dates to compute standard errors. For each metric and node size, we obtain the minimum density of released particles that optimize robustness and computational cost of both metrics when correlation coefficients are above a certain threshold.

The same set of simulations are also used to evaluate how the connectivity metrics are affected by the choice of node size (i.e. the spatial discretization to construct the connectivity matrix). We investigate differences in the spatial patterns of the four metrics when the node size of the network is changed. In fact, for such analysis, Pearson correlation can not be used because it assumes the same number of data-points for comparison, which means it can only be used for pairs of experiments with equal node-sizes.

## Biological parameters

A common use of dispersal models in marine ecology is to evidence the temporal variability of connectivity of a given site (e.g. Qian *et al.* 2015) or to discriminate geographical locations according to their mean connectivity (e.g. Mbaye *et al.* 2015). As such, the starting-date sensitivity study aims at quantifying how connectivity metrics change due to uncertainties in the day of one single spawning event. For this purpose, we compute the Pearson correlation coefficients for each connectivity metric node-by-node comparing pairs of experiments that only differ by their starting dates, maintaining a fixed time-lag (i.e. 2 events separated by a fixed number of days). To do this, we used an ensemble of simulations conducted with particles released every day from 1 September 2010 to 30 November 2010 (total of 91 spawning dates) using a PLD of 30 and 60 days, node size
of 1/2° and 1/8° and the required particle density of 0.518 km⁻² (see table 1). Over the total 91 experiments for a given node size and PLD, we compute average Pearson correlation and its standard error for 51 pairs of starting dates separated by time-lags spanning 1 to 40 days. We expect a decrease in correlation as the time lag between the two dates being compared increases.

Next we explored the sensitivity to spawning frequency. We consider as the reference simulation an hypothetical organism that would spawn every day during 3 months. We choose a daily spawning because is the temporal resolution of the hydrodynamic model and a duration of 3 months to explore a wide range of frequencies. Our objectives are then to evaluate the errors made on the estimation of both the mean and the standard deviation of our four connectivity metric using lower spawning frequencies. The simulations are the same as for the spawning-date sensitivity study (see table 1): connectivity metrics are computed every day from 1 September 2010 to 30 November 2010 (91 spawning dates). As in the previous case, node-size and PLD are fixed when we calculate the temporal averages of the four metrics locally (node-by-node). For all node sizes, PLDs and spawning frequencies we computed the mean relative error (MRE) averaged over all the nodes with respect to the reference value obtained with daily spawning:

$$MRE(f_S) = \sqrt{\frac{1}{N} \sum_{i=1}^{N} \left( \frac{CM_{f_S}^i - CM_{f_D}^i}{CM_{f_D}^i} \right)^2} \quad (1)$$

Where CM is one of the four connectivity metrics, f_D and f_S are the daily frequency and the spawning frequency being tested, respectively. The sum runs over the total N nodes (or geographical sites).

We now explore how our ability to significantly differentiate sites according to their simulated connectivities is affected by the spawning period. We used the same simulations as for the previous sensitivity studies and we fixed PLD, node size and the frequency of spawning to compare the connectivity of all pairs of nodes using the Kruskal-Wallis test (KW). This is a non-parametric version of one-way ANOVA, which tests whether both samples originate from the same distribution (in other words, if the connectivity properties of both nodes can be considered equal). Our reference is the total number of discriminated sites based on an hypothetical organism that spawn every day over 3 months. We then study how the total number of discriminated sites (i.e. the number of pairs of sites which have different connectivity properties according to the KW test) decreases when the spawning period increases for different significance levels.

Finally we test how uncertainties on the PLD would impact the modelled connectivity estimates when considering a single-day spawning event. To do this, we performed simulations releasing particles on 1 September and 1 October every year over 2001-2010 (total of 20 starting dates) with PLD ranging from 5 to 65 days. Tracking times over 65
days are not tested because the 2-dimensional assumptions becomes inaccurate at longer time-scales (Rossi et al. 2014). Particle density was fixed to 0.518 $km^{-2}$ and node-size 1/2° (see table 1). Pearson correlations were computed for the four connectivity metrics by comparing all pairs of experiments differing by their PLDs (but using the same starting dates) and then by averaging the Pearson coefficients. In this way we measure the influence of all the possible positive and negative PLD deviations within our range of PLDs (60 values, and 3600 pairwise comparisons).

## 329 Results

### 330 Numerical parameters

We first investigate the cloud of points on the scatter-plots comparing all nodes of two simulations generated with different densities of particle released (Fig. 2), we observe that most nodes with the highest land ratios depart from the linear relationship. Those nodes are indeed coinciding with the boundaries of the hydrodynamical model (see Fig. 1). Since land ratio is related to initial number of particles per node<fnref n="3" />, outliers are the few marginal nodes which have initially less particles than fully-oceanic nodes. Another issue is related to highly divergent oceanic areas: after a long integration time, the nodes located in such regions may end up with no particle at all (See Supplementary Material, Fig. SM-1). This also produces a sampling error in our metrics. As already mentioned by Dubois et al. (2016), this effect concerns essentially the Alboran Sea and, to a lesser extent, the Ierapetra wind-driven semi-permanent gyre (south of Cyprus).

After further exploratory tests (see also Fig. SM-2 and Fig. SM-3) and to avoid the influence of these issues on our global analyses, we decided in the following to restrict our domain (see the red rectangle in Fig. 1) by excluding nodes whose longitude is < 1°W (removing the Alboran Sea) and > 22°E (removing the Ierapetra gyre), and excluding also the nodes with land ratio greater than 0.8. Note that this exclusion still keeps in the analysis the majority of the coastal nodes.

Sensitivity of the four metrics with respect to density of released particles is shown in figure 3. Instead of displaying Pearson correlation $r$ versus particle density, we plot $1-r$ versus number of particles per node. Thus, decreasing curves indicate that correlation improves with larger particle density. For all connectivity metrics, increasing the number of particles released per node returns higher $r$, so that $1-r$ goes to 0. Moreover, the scaling of $(1-r)$ with the number of particles released per node for the connectivity metrics LR, SR and $SS_s$ is approximately a power law with exponent -1 independently

<footnote n="3">³Note that the seeding of particles only occurs in the vicinity of all effective grid points of the velocity field. In this way, we do not seed particles that will not be advected during the simulation (to avoid over-estimating retention).</footnote>
<table><caption>Table 2: Minimum particle density and number of particles per node required to obtain a robust result for $SS_d$ choosing a threshold of 0.95 in $r$. The thresholding procedure for the other metrics ($LR$, $SR$ and $SS_s$) returns the single value of ≈ 100 particles per node, for any node-size.</caption><thead><tr><th colspan="3">Threshold (r ≥ 0.95) for SS<sub>d</sub></th></tr><tr><th>Node size (°)</th><th>Particle density (km<sup>-2</sup>)</th><th>Number of particles per node</th></tr></thead><tbody><tr><td>1</td><td>0.518</td><td>6400</td></tr><tr><td>3/4</td><td>0.518</td><td>3600</td></tr><tr><td>1/2</td><td>0.518</td><td>1600</td></tr><tr><td>1/4</td><td>1.037</td><td>800</td></tr><tr><td>1/8</td><td>1.037</td><td>200</td></tr><tr><td>1/10</td><td>1.556</td><td>≈ 190</td></tr><tr><td>1/12</td><td>1.556</td><td>≈ 130</td></tr></tbody></table>

of node-size (black line in Fig. 3). However, in the case of $SS_d$ this dependence is more complex: it still decreases but slowly than the other metrics.

Based on this figure, one can estimate the minimum particle density required to obtain a robust result for each connectivity metric. We used a threshold of 0.95 for $r$ ($1-r = 0.05$ in Fig. 3) so that we consider the connectivity estimates to be insensitive to the number of particles released when $r$ is greater than this value. We found a threshold of ≈ 100 particles per node for all node sizes when considering the metrics that are computed from the connectivity matrix ($LR$, $SR$ and $SS_s$). Concerning $SS_d$, the thresholding procedure shows a dependence on the node size: values are reported in table 2.

Now we investigate the spatial variability looking at maps of the four metrics computed at a given time while varying the node size. A visual comparison between panels of figures 4, 5 and 6 is informative. Spatial patterns are finer for the high resolution discretization but large scale patterns still remain in lower resolutions. However, as the node size increased, the spatial variability decreased in exchange metrics ($SS_s$ and $SS_d$), on the contrary, spatial variability decreases for higher grid resolutions in retention metrics ($LR$ and $SR$) since the number of nodes with zero value increases.

## Biological parameters

Results of spawning date sensitivity study for node sizes of 1/2° and 1/8°, using the two PLDs 30 and 60 days are shown in figure 7. Overall, Pearson correlation coefficients decrease when the time-lag between starting dates increases. For all metrics, one can distinguish two regimes with a quite rapid decrease during the first 8 to 15 days and then a rather slow linear decrease of the correlation until time lag matches with PLD value. At time lags higher than PLD the correlation stabilizes, however, our time-lag range only allows us to observe this for the PLD of 30 days. Note also that while the retention indices ($LR$ and $SR$) remain quite robust (coefficients spanning 0.55 - 1), the exchange
indices reveal higher dependence ($SS_s$ and $SS_d$) with respect to spawning times with $r$ falling down to 0.2 or lower for lag-times of 40 days. There are not large differences between both PLDs but it is worth noting that shorter PLDs are more sensitive to the spawning date for all metrics. We observe also that for larger node sizes the robustness to spawning-date variations increases for exchange metrics but decreases for retention ones.

Results of the spawning-frequency sensitivity study are shown in figure 8 and 9. Overall, the errors made when estimating the mean connectivity and its fluctuation as compared with the reference case of daily spawning increase with spawning period (Fig. 8). Mean Relative Error (MRE) is always larger (about one order of magnitude) in retention metrics than exchange ones. PLD has very little influence on the MRE for all the frequencies examined, but larger node sizes tend to produce lower mean relative errors. For instance, if one uses a node size of 1/8° and a spawning period of 10 days over 3 months to study our hypothetical species, the mean connectivity estimates of retention will have errors (with respect to the reference daily spawning) of about 5% while exchange metrics will have errors of 50 to 100%, independently of its PLD.

When investigating the percentage of pairs of sites that can be statistically distinguished by their connectivity-metric distributions we see that it decreases with the spawning periodicity (Fig. 9 plots this percentage relative to the value for daily spawning). Note that this percentage decreases substantially slower for the retention metrics than for the exchange metrics, while $SS_d$ has the faster decay. In this case, the PLD has an influence: long PLDs resulted in an increase (marked in $SS_s$ and $SS_d$) of the number of discriminated pairs (see Fig. SM-6).

Results of the PLD sensitivity study are shown in figure 10. The four connectivity metrics reveal different behaviour. The exchange proxies $SS_s$ and $SS_d$ are considerably less robust to variations of the PLD than the metrics measuring retention (see the area enclosed between the $r = 0.9$ black contours, Fig. 10). Another interesting result is that uncertainties of a few days for long PLDs have less impact than a few-days uncertainty for small PLDs, as evidenced by the conical structures visible of isolines in all subplots of Fig. 10, especially marked within LR and SR.

# Discussion and conclusions

## Numerical parameters

We analysed the impact of the number of released particles on our different connectivity metrics and found two different behaviours: $SS_d$ improves with the number of particles slower than the other metrics. Moreover its scaling (i.e. the way in which the correlation with the reference case increases for larger particle number) depends also on the node-size, in contrast with the scaling of LR, SR and $SS_s$ which depends only on the initial
number of particles per node. Within our setup we find that 100 particles per node are sufficient to ensure the robustness of *LR*, *SR* and *SS*<sub>s</sub> independently of node size. It is worth noting that this precise number may depend on other factors not included in this study: the characteristics of the underlying velocity field or the implementation of larval behavior (e.g. active vertical migration). However, the robustness threshold for *SS*<sub>d</sub> depends on node side and remains greater than 100 particles (table 2) due to its different scaling. Simons *et al*. (2013) also found a minimum number of particles which guarantees robustness in their model predictions. They found an exponential scaling of robustness with number of particles released but note that their results are not completely comparable to ours due to methodological differences. Our analyses focused on a PLD 60 days, however we stress that for shorter PLDs the recommended number of particles per node would be lower than the threshold we propose. This is because the stochasticity of trajectories and the complexity of the circulation patterns experienced by particles decreased with tracking time considering the passive tracer assumption. Moreover, the consistent threshold of 100 particles per node we found is independent on node size. The exact threshold for number of particles to be released would essentially depend on the characteristics of the flow field (and corresponding ocean model) but the general scaling laws we found (functional form of the improvement of correlation with particle number) still hold for any off-line particle model.

Maps of the metric fields showed that more detailed patterns can be revealed using small nodes for both retention and exchange. Concerning *LR* and *SR* (Fig. 4), smaller nodes lead to non-zero values (red colors) essentially in nearshore and shelf regions, while those same retentive regions are less intense when considering larger nodes (i.e. which extend further away from continental shelves, thus encompassing the highly dispersive character of energetic boundary currents), consequently spatial variability increases with larger node sizes. This suggests that to obtain reliable estimates of larval retention, the scales of discretization (node-size) should be at least of the order of, or lower, than the mean width of the continental shelves. This is consistent with Brochier *et al.* (2011); Mbaye *et al.* (2015) who uses high resolution computation to evidence high retention rates nearshore and over the continental shelves. Regarding the proxies of exchange, we stress that in large nodes exchange processes result to be relatively less important than in small nodes, therefore variability of values decrease. Note however that the relative range of values and the spatial patterns are similar among different node-sizes. The high robustness of *SS* patterns when increasing the coarse-graining of the system is related to the robustness of Lyapunov exponents (Hernandez-Carrasco *et al.* 2011). Indeed, Ser-Giacomi *et al.* (2015a) demonstrated a relationship among such exponents and network measures used to compute the *SS*<sub>d</sub> metrics. Note however that smaller nodes need higher computation costs, especially when interested in *SS*<sub>d</sub> which return the highest threshold of density particles per area (Table 2).
The choice of the node size has biological implications. As seen previously, it must be determined according to scales (e.g. basin versus local) and processes (retention or exchange) related to the biological hypothesis to be tested. For instance, different parameter choices must be made if studying the dispersal of juvenile sea turtles across the entire Atlantic (Putman & He 2013) or among a smaller region like the Balearic archipelago. Secondly, one must consider the typical length-scales of fragmentation of the preferred habitat of the target species. For instance, small node sizes permit a better consideration of highly fragmented habitats (e.g. benthic habitats) and improved detection of small-scale dispersal patterns. However we also found that exchange metrics are highly dependent on the node-size, which affects the quantification of import/export from neighbouring patches. The importance of the choice of length-scales becomes even more relevant in seascape ecology (i.e. where no clear boundaries delimitate habitats). Indeed, Alvarez-Berastegui et al. (2014) demonstrate that a combination of hydrographic variables and their spatial gradients improve their ability to model pelagic fish distribution only if they are calculated at the appropriate spatial scales. This means that to incorporate the quality of pelagic habitats within modelled connectivity estimates, the resolution of hydrographic variables and the spatial-scales at which they are averaged should be determined accordingly to the node size.

These conclusions have broad implications as they will help any study using particle-tracking model to choose wisely the density of initial particles and the size of discretization in accord with the objectives of the study and the range of parameters considered.

## Biological parameters

The sensitivity study to spawning date highlighted that retention metrics are always more robust than exchange ones. We also found two regimes for all PLDs and node-sizes studied. This can be explained by temporal correlations of Lagrangian velocities. In the first regime, velocities are correlated among the reference simulation and lagged ones, this generates a fast decay of correlations, since it produces oriented changes in spatial patterns. The second regime is characterized by uncorrelated velocities which causes the slow down of pairwise correlations due to random variations of spatial patterns. The time lag at the transition between both regimes can be estimated by the Lagrangian correlation time (Taylor 1922), i.e. the time at which the Lagrangian velocities start to be uncorrelated. We found here that the second regime starts approximately at time-lag of 10 days, which is greater than the Lagrangian correlation times of about 2-5 days in the ocean (Poulain & Zambianchi 2007; Berti et al. 2011). This discrepancy might be due to that fact that for time lags larger than the Lagrangian correlation time, velocities can still be correlated at particular locations, especially due to mesoscale features (e.g. larvae trapped within eddies) (Elhmaidi et al. 2010; Brochier et al. 2011). Looking at
PLDs of 30 days, the existence of a third regime is suggested when the time-lag becomes greater than the PLD (correlation coefficients are constant and low). It could mean that the pair of experiments are statistically equivalent, therefore there is no causality between both simulations that do not overlap temporally. We only observe this for one tested PLD because the maximum time-lag used here is kept below 60 days. Therefore additional studies, preferentially in three-dimensions, are needed to further explore longer time-scales.

Correlation of connectivity metrics obtained with lagged spawning dates revealed that they are significantly more robust for longer PLD. Therefore, the uncertainties about the spawning dates for a short-PLD species will largely affect the results, whereas the same level of uncertainty for a long-PLD species will have limited impacts.

The influence of node size on the robustness of spawning date depends on the metric considered. Exchange metrics ($SS_s$ and $SS_d$) are more robust when node-size increases, independently of the time-lag. On the contrary, retention metrics ($LR$ and $SR$) are more robust when node size decreases. As we discussed previously (sect. *Numerical parameters*), node size rules the spatial variability of retention and exchange: larger nodes return more variable retention processes, hence retention metrics are less robust, and smaller nodes increment the rate of exchanges resulting in a weaker robustness for exchange metrics.

The sensitivity to spawning frequency was studied from two perspectives: we first studied how decreasing the spawning frequency with respect a reference case of daily spawning increases the errors in the estimates of connectivity proxies for each node; then we investigated how the ability to differentiate distinct geographical sites by their connectivity properties decreases when spawning period is increased. We also found that, with increasing spawning periodicity, errors grow slower in retention proxies than in exchange proxies (about one order of magnitude). This is consistent with previous results and it suggest that changes in retention processes are primarily spatially-driven, i.e. controlled by the bathymetry that influences the local currents. In contrast, temporal variability is very high in exchange metrics, so that high frequency must be considered to properly describe them. We also found that PLD has very little influence on metrics errors when the spawning period is lower than the PLD itself. Furthermore, increasing node size always decrease errors in both retention and exchange metrics.

We showed that increase in the spawning period affects negatively the statistical discrimination of sites with distinct connectivity properties. This is in line with Kough & Paris (2015) who found that spawning periodicity exerts a control on the number and persistence of network connections among distant sub-populations. This negative influence is considerably less intense in retention metrics than exchange ones. It reinforces our previous statement: any retention proxies display lower temporal variability than exchange metrics. By combining our conclusions derived from Fig. 8 and 9, we found that,
when increasing the periodicity, small relative errors in retention metrics have relatively large negative impacts on the discriminations of sites. In contrast, large relative errors made on exchange metrics could still allow to discriminate sites properly. This has to do with the distinct levels of spatial and temporal variability of retention and exchange metrics (Qian et al. 2015), which are directly related to the intrinsic properties of the velocity field, and thus of the hydrodynamical model.

This spawning-frequency study has biological implications: if a given species is known to spawn continuously over a 2-months period, similar mean connectivity estimates will be obtained when considering daily spawning or periods of 5 to 10 days. In other words, for a species known to spawn repeatedly over long period of a few months, using a periodicity of 5 to 10 days will still allow discriminating sites by their mean rates of larval retention and exchange. This is consistent with Huggett et al. (2003) who reported similar transport success for a PLD span over 15 days to 2 months when using periods ranging from 1 to 10 days. Snyder et al. (2014) argued that connectivity fluctuations (i.e. the temporal variability of ocean currents echoing on the variability of connectivity) would affect more the seasonal spawners than the year-round spawners. Our results do not support this hypothesis since the 10 days threshold would apply to most oceanic regions and seasons. Nevertheless, it is also clear that the level of connectivity fluctuations highly depends on the locations and periods considered (e.g. Fig. SM-4, SM-5). Indeed, there exist oceanic areas (e.g. the lagoon system of Cuif et al. 2014) for which higher spawning frequencies (i.e. spawning periods smaller than one day) are required to take into account rapidly changing flow fields (e.g. tidal currents, sea-breeze induced circulation, short-living sub-mesoscale fronts and eddies, etc...). Besides, our study suggests that the most precise connectivity estimates for given pairs of location/period could be attained by assessing *a-priori* the characteristic spatial and temporal variability scales of the flow and by selecting parameters accordingly. This is because, despite the stochasticity of the ocean circulation, there are locations and periods which have relatively stable dynamical properties, such as its level of turbulence (that could be evaluated through its eddy kinetic energy or Lyapunov fields) and its divergent/convergent character, thus resulting in quite constant connectivity patterns (see Dubois et al. (2016); see also Fig. SM-4a, SM-5a for examples for constant connectivity metrics and SM-4b and SM-5b for examples of highly variable ones).

Concerning the sensitivity to the PLD, we observe that correlations decrease when the PLD is shifted from a reference PLD. However, the rate of decrease for retention metrics revealed stronger dependence on the reference PLD than for exchange diagnostics. While it is well known that long PLDs increase mean dispersal distances of propagules and decrease their retention (Garavelli et al. 2014) (Fig. SM-5a), the effects on exchange diagnostics may be not straightforward (Fig. SM-5b). We also found that relative errors in connectivity estimates due to uncertainties on the PLD are smaller for long PLD
(Fig. 10). This means that if the mean PLD of a given species is 50 days ±10 days, the computation of retention indices will be very similar using 60 or 40 days tracking time (the corresponding pair or *r*-values falls within the black contours of Fig. 10). In contrast, for a given species with a mean PLD of 15 days ±10 days, retention rates will be significantly different for PLDs of 5 or 25 days (outside the black contours in Fig. 10). This is consistent with Simons *et al.* (2013) who found that the sensitivity of particle density distribution to PLD is much larger for short tracking times of 15 days or less than for longer tracking times of 20 days or more. We also demonstrated that this robustness effect does not hold for exchange metrics. Note however that for PLD uncertainties lower than the characteristic time-scale of the flow (i.e. ± a few days around the reference PLD), exchange metrics remain constant because the convergence of the flow is almost unchanged at these time-scales.

## Other factors not considered

While some species-specific factors were not included in our sensitivity studies, they have the potential to affect connectivity estimates and their robustness. We did not explicitly test the sensitivity of our metrics to the depth of dispersion because it would depend on the location and season of the oceanic region of interest, as well as on the mixing parametrization of the hydrodynamical model. In the Baltic sea, Corell *et al.* (2012) found that the depth distribution of larvae critically affects their dispersal. Along the Chilean coasts Garavelli *et al.* (2014) found that retention increases with depth since currents are generally weaker and often shoreward. Simons *et al.* (2013) suggested that the robustness of the depth of dispersion on connectivity estimates is proportional to the strength of vertical mixing. In other words, summer well-stratified waters are usually dominated by baroclinic currents so that connectivity estimates are sensitive to the depth of dispersion. In contrast, wintertime deep-mixed layers waters and weak stratification favour more barotropic flow in which connectivity estimates are less sensitive to the depths of dispersion.

In the Mediterranean Sea, Torres *et al.* (2011) observed in spring that most fish larvae are concentrated in the surface mixed-layer, that is between 20 and 50 m. This is consistent with d'Ortenzio *et al.* (2005) who reported surface mixed-layer of less than 50 m in summer and about 100 m during winter or even more in specific areas of deep-water formation. It suggests that connectivity estimates during summer might be very similar in the top 50 m, extending down to 100 m in winter. This is consistent with Mariani *et al.* (2010) who argued that mesoscale processes (whose magnitude is governed by atmospheric seasonal forcing) control aggregation and dispersion processes in the north-west Mediterranean sea. Note that from a more global perspective, Sundby & Kristiansen (2015) proposed a simplified classification of fish eggs as pelagic (0-50 m), mesopelagic
(150-250 m) and bottom eggs (250-350 m, or the bottom of the shelf). Since dispersal at various depths can be examined by selecting the corresponding vertical layer of the model, it will become possible to test this effect globally.

Another factor that was not tested here but that might affect connectivity estimates is the representation of small-scale (sub-mesoscale and below) structures in the simulated velocity fields. Indeed, Putman & He (2013) indicated that the use of non-eddy-resolving model outputs to compute drifting trajectories yields predictions that are inconsistent with observations. To circumvent this issue, it is common to add a random component (or random walk) to the particle trajectories to reproduce this non-resolved turbulence. Qin et al. (2014) found that adding a random-walk term to the trajectories calculated with an eddy-resolving model does not mitigate the errors introduced by the temporal averaging, at least for regional to large-scale focus. By adding noise to trajectories or by smoothing a given velocity field, Hernandez-Carrasco et al. (2011) investigate the effect of the effective horizontal resolution of the simulated currents and found small impacts on the computation of Lyapunov exponents. We expect the same levels of robustness for our connectivity metrics because network measures (such as degrees, which allows the computation of $SS_d$) are physically related to the Lyapunov exponents field (Ser-Giacomi et al. 2015a). Concerning the vertical dimension, Simons et al. (2013) found that adding vertical sub-grid scale mixing does not alleviate the sensitivity of trajectories to the depth of dispersion. Although the two latter studies advocate for a small impacts of sub-grid mixing to connectivity estimates, further studies are needed, especially considering the growing development of very high-resolution hydrodynamical models.

Last, the implementations of complex larval behaviour, mortality or success of settlement (e.g. Lett et al. 2008; Paris et al. 2013) as well as three-dimensional trajectories (Kool & Nichol 2015) are not studied here but are envisaged for future development. While Garavelli et al. (2016) reported that vertical migration is a biological factor that significantly influences connectivity patterns, it is likely it would be so primarily in systems characterized by strong vertical shear (e.g. coastal upwelling Queiroga et al. 2007; Brochier et al. 2011). Another relevant process is the propagule buoyancy (e.g. changing mass density of eggs during its early-life stages) which has been found to significantly affect transport success of anchovy eggs (Parada et al. 2003) but not in the larvae of a sessile gastropod Garavelli et al. (2012). This example clearly suggests that the effect of those fine biological parametrization on connectivity are highly species-specific.

## General conclusions

We have analysed the global robustness and sensitivity of four connectivity metrics measuring the exchange and retention of larvae derived from the LFN framework. The most important numerical and bio-physical parameters were extensively tested over large
ranges. Despite the fact that our proxies are sensitive at different degrees to all parame-
ters, we evidenced general relationships that demonstrated the robustness and identified
the limitations of such particle-tracking models. In particular, we found an optimum par-
ticle density that ensures the best ratio computational cost/information obtained. We
assessed the impacts of the scale of discretization (node size) and we provided hints on
how to select them adequately according to the study objectives. The consequences of un-
certainties about the PLD and the modality of spawning were systematically explored and
related to the spatio-temporal variability of the underlying velocity field. Such variability
can be characterized by the Lagrangian correlation time: a measure of the time-scale at
which fluid velocities become uncorrelated. We found retention metrics are always more
robust to the temporal variability of ocean currents than exchange metrics.

Our results have implications: they help designing suitable experiments with the LFNs
framework to address specific ecological questions and allow evaluating the reliability one
may give to these model estimates. These results, the methodologies we employed and
more generally the LFN modelling framework are equally applicable to the study of
other passively-transported materials (oil spills, micro-plastics,...) from a Lagrangian
perspective. Future perspectives are to relate these robustness results with physical
quantities in order to extend our results to different regions or conditions and improve
the results including biological traits and behavior.

References

Alvarez-Berastegui, D., Ciannelli, L., Aparicio-Gonzalez, A., Reglero, P., & Hidalgo, M. et al. 2014. Spatial Scales, Means and Gradients of Hydrographic Variables Define Pelagic Seascapes of Bluefin and Bullet Tuna Spawning Distribution. *PLoS ONE*, **9**(10), e109338.

Andrello, M., Mouillot, D., Beuvier, J., Albouy, C., Thuiller, W., & Manel, S. 2013. Low connectivity between Mediterranean marine protected areas: a biophysical modeling approach for the dusky grouper *Epinephelus marginatus*. *PLoS ONE*, **8**(7), e68564.

Berti, S., Alves Dos Santos, F., Lacorata, G., & Vulpiani, A. 2011. Lagrangian Drifter Dispersion in the Southwestern Atlantic Ocean. *Journal of Physical Oceanography*, **41**, 1659–1672.

Brochier, T., Lett, C., & Fréon, P. 2011. Investigating the ‘northern Humboldt para-
dox’ from model comparisons of small pelagic fish reproductive strategies in eastern
boundary upwelling ecosystems. *Fish and Fisheries*, **12**(1), 94–109.

Calò, A., Félix-Hackradt, F.C., Garcia, J., Hackradt, C.W., Rocklin, D., Treviño Otón, J.,
& García-Charton, J.A. 2013. A review of methods to assess connectivity and dispersal
679

680

between fish populations in the Mediterranean Sea. *Advances in Oceanography and Limnology*, **4**(2), 150-175.

681

682

683

Capet, X., McWilliams, J-C., Molemaker, M.J., & F., Shchepetkin. 2008. Mesoscale to Submesoscale Transition in the California Current System. Part I: Flow Structure, Eddy Flux, and Observational Tests. *Journal of Physical Oceanography*, **38**, 29-43.

684

685

686

Carr, M.H., Neigel, J.E., Estes, J.A., Andelman, S., Warner, R.R., & Largier, J.L. 2003. Comparing marine and terrestrial ecosystems: implications for the design of coastal marine reserves. *Ecological Applications*, **13**(1), 90-107.

687

688

689

Clayton, S., Dutkiewicz, S., Jahn, O., & Follows, M.J. 2013. Dispersal, eddies, and the diversity of marine phytoplankton. *Limnology and Oceanography: Fluids and Environments*, **3**, 182-197.

690

691

692

Corell, H., Moksnes, P.O., Engqvist, A., Doos, K., & Jonsson, P. 2012. Depth distribution of larvae critically affects their dispersal and the efficiency of marine protected areas. *Marine Ecology Progress Series*, **467**, 29-46.

693

694

Cowen, R.K., & Sponaugle, S. 2009. Larval dispersal and marine population connectivity. *Annual Review of Marine Science*, **1**, 443-466.

695

696

697

698

Cuif, M., Kaplan, D.M., Lefèvre, J., Faure, V.M., Caillaud, M., Verley, P., Vigliola, L., & Lett, C. 2014. Wind-induced variability in larval retention in a coral reef system: A biophysical modelling study in the South-West Lagoon of New Caledonia. *Progress in Oceanography*, **122**, 105-115.

699

700

d'Ortenzio, F., & d'Alcalá, M. R. 2009. On the trophic regimes of the Mediterranean Sea: A satellite analysis. *Biogeosciences*, **6**, 1-10.

701

702

703

704

d'Ortenzio, F., Iudicone, D., de Boyer Montegut, C., Testor, P., Antoine, D., Marullo, S., Santoleri, R., & Madec, G. 2005. Seasonal variability of the mixed layer depth in the Mediterranean Sea as derived from in situ profiles. *Geophysical Research Letters*, **32**, L12605.

705

706

707

708

Dubois, M., Rossi, V., Ser-Giacomi, E., Arnaud-Haond, S., López, C., & Hernández-García, E. 2016. Linking basin-scale connectivity, oceanography and population dynamics for the conservation and management of marine ecosystems. *Global Ecology and Biogeography*, **25**(5), 503-515.

709

710

Duputie, A., & Massol, F. 2013. An empiricist's guide to theoretical predictions on the evolution of dispersal. *Interface Focus*, **3**, 20130028.
711

712

Elhmaidi, D., Nefzi, H., Carton, X., & Lili, T. 2010. Particle Dispersion in the Western Mediterranean Basin. *The Open Oceanography Journal*, **4**, 137–143.

713

714

Froyland, G., Padberg, K., England, M.H., & Treguier, A.-M. 2007. Detection of coherent oceanic structures via transfer operators. *Physical Review Letters*, **98**(22), 224503.

715

716

717

Froyland, G., Horenkamp, C., Rossi, V., Santitissadeekorn, N., & Sen Gupta, A. 2012. Three-dimensional characterization and tracking of an Agulhas Ring. *Ocean Modelling*, **52-53**, 69–75.

718

719

720

Garavelli, L. Grüss, A., Grote, B., Chang, N., Smith, M., Verley, P., Stenevik, E.K., Kaplan, D.M., & Lett, C. 2012. Modeling the Dispersal of Cape Hake Ichthyoplankton. *Journal of Plankton Research*, **34**, 655–669.

721

722

723

724

Garavelli, L., Kaplan, DM., Colas, F., Stotz, W., Yannicelli, B., & Lett, C. 2014. Identifying appropriate spatial scales for marine conservation and management using a larval dispersal model: The case of Concholepas concholepas (loco) in Chile. *Progress in Oceanography*, **124**, 42–53.

725

726

727

Garavelli, L., Colas, F., Verley, P., Kaplan, D.M., Yannicelli, B., & Lett, C. 2016. Influence of Biological Factors on Connectivity Patterns for Concholepas concholepas (loco) in Chile. *PloS ONE*, **11**, 1–22.

728

729

730

Guizien, K., Belharet, M., Marsaleix, P., & Guarini, J.M. 2012. Using larval dispersal simulations for marine protected area design: Application to the Gulf of Lions (north-west Mediterranean). *Limnology and Oceanography*, **57**(4), 1099.

731

732

733

Harrison, C.S., Siegel, D.A., & Mitarai, S. 2013. Filamentation and eddy–eddy interactions in marine larval accumulation and transport. *Marine Ecology Progress Series*, **472**, 27–44.

734

735

736

Hernandez-Carrasco, I., López, C., Hernández-García, E., & Turiel, A. 2011. How reliable are Finite-Size Lyapunov Exponents for the assessment of ocean dynamics? *Ocean Modelling*, **36**, 208–218.

737

738

739

Huggett, J., Freon, P., Mullon, C., & Penven, P. 2003. Modelling the transport success of anchovy Engraulis encrasicolus eggs and larvae in the southern Benguela: The effect of spatio-temporal spawning patterns. *Marine Ecology-Progress Series*, **250**, 247–262.

740

741

742

743

Huret, M., Petitgas, P., & Woillez, M. 2010. Dispersal kernels and their drivers captured with a hydrodynamic model and spatial indices: A case study on anchovy (Engraulis encrasicolus) early life stages in the Bay of Biscay. *Progress In Oceanography*, **87**(1-4), 6–17.
744
745

Jones, G.P., Srinivasan, M., & Almany, G.R. 2007. Population connectivity and conservation of marine biodiversity. *Oceanography*, **20**, 42–53.

746
747

Jönsson, B.F., & Watson, J.R. 2016. The timescales of global surface-ocean connectivity. *Nature Communications*, **7**, 1–6.

748
749
750

Kool, J.T., & Nichol, S.L. 2015. Four-dimensional connectivity modelling with application to Australia's north and northwest marine environments. *Environmental Modelling & Software*, **65**, 67–78.

751
752

Kool, J.T., Moilanen, A., & Treml, E.A. 2013. Population connectivity: recent advances and new perspectives. *Landscape Ecology*, **28**, 165–185.

753
754

Kough, A. S., & Paris, C. B. 2015. The influence of spawning periodicity on population connectivity. *Coral Reefs*, **34**(3), 753–757.

755
756
757

Lett, C., Verley, P., Mullon, C., Parada, C., Brochier, T., Penven, P., & Blanke, B. 2008. A Lagrangian tool for modelling ichthyoplankton dynamics. *Environmental Modelling & Software*, **23**(9), 1210–1214.

758
759

Lowe, W.H., & Allendorf, F.W. 2010. What can genetics tell us about population connectivity? *Molecular Ecology*, **19**, 3038–3051.

760
761
762

MacPherson, E., & Raventós, N. 2006. Relationship between pelagic larval duration and geographic distribution in Mediterranean littoral fishes. *Marine Ecology Progress Series*, **327**, 257–265.

763

Mahadevan, A. 2014. Eddy effects on biogeochemistry. *Nature*, **506**, 168–169.

764
765
766

Mancho, A.M., Hernandez-Garcia, E., Small, D., Wiggins, S., & Fernandez, V. 2008. Lagrangian transport through an ocean front in the North-Western Mediterranean Sea. *Journal of Physical Oceanography*, **38**, 1222–1237.

767
768
769

Mariani, P., MacKenzie, B.R., Iudicone, D., & Bozec, A. 2010. Modelling retention and dispersion mechanisms of bluefin tuna eggs and larvae in the northwest Mediterranean Sea. *Progress in Oceanography*, **86**(1), 45–58.

770
771
772

Mbaye, B.C., Brochier, T., Echevin, V., Lazar, A., Lévy, M., Mason, E., Gaye, A.T., & Machu, E. 2015. Do Sardinella aurita spawning seasons match local retention patterns in the Senegalese–Mauritanian upwelling region? *Fisheries Oceanography*, **24**, 69–89.

773
774
775

Mitarai, S., Siegel, D.A., Watson, J.R., Dong, C., & McWilliams, J.C. 2009. Quantifying connectivity in the coastal ocean with application to the Southern California Bight. *Journal of Geophysical Research - Oceans*, **114**(C10026).
776

777

778

Oddo, P., Adani, M., Pinardi, N., Fratianni, C., Tonani, M., & Pettenuzzo, D. 2009. A nested Atlantic-Mediterranean Sea general circulation model for operational forecasting. *Ocean Science Discussions*, **6**(2), 1093–1127.

779

780

Palumbi, S.R. 2003. Population genetics, demographic connectivity, and the design of marine reserves. *Ecological Applications*, **13**(sp1), 146–158.

781

782

783

784

Parada, C., van der Lingen, C.D., Mullon, C., & Penven, P. 2003. Modelling the effect of buoyancy on the transport of anchovy (Engraulis capensis) eggs from spawning to nursery grounds in the southern Benguela: an IBM approach. *Fisheries Oceanography*, **12**, 170–184.

785

786

787

Paris, C.B., Chérubin, L.M., & Cowen, R.K. 2007. Surfing, spinning, or diving from reef to reef: effects on population connectivity. *Marine Ecology Progress Series*, **347**, 285–300.

788

789

790

Paris, C.B., Helgers, J., van Sebille, E., & Srinivasan, A. 2013. Connectivity Modeling System: A probabilistic modeling tool for the multi-scale tracking of biotic and abiotic variability in the ocean. *Environmental Modelling and Software*, **42**, 47–54.

791

792

793

Poulain, P.M., & Zambianchi, E. 2007. Surface circulation in the central Mediterranean Sea as deduced from Lagrangian drifters in the 1990s. *Continental Shelf Research*, **27**, 981–1001.

794

795

796

Putman, N.F., & He, R. 2013. Tracking the long-distance dispersal of marine organisms: sensitivity to ocean model resolution. *Journal of the Royal Society Interface*, **10**(20120979).

797

798

Qian, H., Li, Y., He, R., & Eggleston, D. B. 2015. Connectivity in the Intra-American Seas and implications for potential larval transport. *Coral Reefs*, **34**(2), 403–417.

799

800

801

Qin, X., van Sebille, E., & Sen Gupta, A. 2014. Quantification of errors induced by temporal resolution on Lagrangian particles in an eddy-resolving model. *Ocean Modelling*, **76**, 20–30.

802

803

804

805

Queiroga, H., Dos-Santos, A., Dubert, J., Gonzalez-Gordillo, I., Paula, J., Peliz, A., & Santos, A.M.P. 2007. Oceanographic and behavioural processes affecting invertebrate larval dispersal and supply in the western Iberia upwelling ecosystem. *Progress in Oceanography*, **74**, 174–191.

806

807

808

809

Reglero, P., Ciannelli, L., Alvarez-Berastegui, D., Balbín, R., López-Jurado, J.L., & Alemany, F. 2012. Geographically and environmentally driven spawning distributions of tuna species in the western Mediterranean Sea. *Marine Ecology Progress Series*, **463**, 273–284.
810

811

812

Rossi, V., López, C., Sudre, J., Hernández-García, E., & Garçon, V. 2008. Comparative study of mixing and biological activity of the Benguela and Canary upwelling systems. *Geophysical Research Letters*, **35**, L11602.

813

814

815

Rossi, V., Van Sebille, E., Sen Gupta, E., Garçon, V., & England, M.H. 2013. Multi-decadal projections of the surface and interior pathways of the Fukushima Cesium-137 radioactive plume. *Deep Sea-Research I*, **80**, 37–46.

816

817

818

Rossi, V., Ser-Giacomi, E., López, C., & Hernández-García, E. 2014. Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves. *Geophysical Research Letters*, **41**(8), 2883–2891.

819

820

Ser-Giacomi, E., Rossi, V., Lopez, C., & Hernandez-Garcia, E. 2015a. Flow networks: A characterization of geophysical fluid transport. *Chaos*, **25**, 036404.

821

822

823

Ser-Giacomi, E., Vasile, R., Hernández-García, E., & López, C. 2015b. Most probable paths in temporal weighted networks: An application to ocean transport. *Physical Review E*, **92**(1), 012818.

824

825

Shanks, A.L. 2009. Pelagic larval duration and dispersal distance revisited. *The Biological Bulletin*, **216**(3), 373–385.

826

827

828

Simons, R.D., Siegel, D.A., & Brown, K.S. 2013. Model sensitivity and robustness in the estimation of larval transport: A study of particle tracking parameters. *Journal of Marine Systems*, **119–120**, 19–29.

829

830

Snyder, R.E., Paris, C.B., & Vaz, A.C. 2014. How much do marine connectivity fluctuations matter? *The American Naturalist*, **184**, 523–530.

831

832

Sundby, S., & Kristiansen, T. 2015. The Principles of Buoyancy in Marine Fish Eggs and Their Vertical Distributions across the World Oceans. *PLoS ONE*, **10**(10), e0138821.

833

834

Taylor, G.I. 1922. Diffusion continuos mouvements. *Proceedings of The London Mathematical Society*, **20**, 196–212.

835

836

837

Torres, A.P., Reglero, P., Balbin, R., Urtizberea, A., & Alemany, F. 2011. Coexistence of larvae of tuna species and other fish in the surface mixed layer in the NW Mediterranean. *Journal of Plankton Research*, **33**(12), 1793–1812.

838

839

840

841

Treml, E.A., Roberts, J.J., Chao, Y., Halpin, P.N., Possingham, H.P., & Riginos, C. 2012. Reproductive output and duration of the pelagic larval stage determine seascape-wide connectivity of marine populations. *Integrative and Comparative Biology*, **52**(4), 525–537.
842

843

844

Vaz, A.C., Richards, K.J., Jia, Y., & Paris, C.B. 2013. Mesoscale flow variability and its impact on connectivity for the island of Hawaii. *Geophysical Research Letters*, **40**, 332–337.

845

846

847

Wood, S., Paris, C.B., Ridgwell, A., & Hendy, E.J. 2014. Modelling dispersal and connectivity of broadcast spawning corals at the global scale. *Global Ecology and Biogeography*, **23**(1), 1–11.

848

## Supplementary Material

849

850

The following supplementary material is available at ICES-JMS online: Supplementary figures SM-1 through SM-6.

851

## Acknowledgements

852

853

854

855

856

857

858

859

860

861

862

The authors would like to thank M. Hidalgo for fruitful discussions. The authors acknowledge supports from Ministerio de Economía y Competitividad and Fondo Europeo de Desarrollo Regional through the LAOP project (CTM2015-66407-P, MINECO/FEDER) and through a Juan de la Cierva – Incorporación fellowship (IJCI-2014-22343) granted to V.R., from ECs Marie-Curie ITN program (FP7-PEOPLE-2011-ITN) through the LINC project (289447) and from MISTRALS ENVI-Med through the HYDROGEN-CONNECT project. This work has received support under the program “Investissements d’Avenir” launched by the French Government and implemented by ANR with the references ANR-10-LABX-54 MEMOLIFE and ANR-11-IDEX-0001-02 PSL Research University. The simulated velocity field used in these studies was generated by Copernicus (www.copernicus.ocean.eu).
<figure><img src="image_3.png" /><figcaption>Figure 2: Scatter-plots for each connectivity metrics ( (a) Local-Retention, (b) Self-Recruitment, (c) Source-Sink strength and (d) Source-Sink degree) among two LFN which only differ in particle density (0.518 and 5.176 km⁻²). The other constant parameters are node-size=1/8°, depth=12 m, starting date is 16 October 2010 and PLD=60 days. Colour represents the value of land ratio in the corresponding node.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Initial particle density sensitivity tests for each connectivity metric: (a) *LR*, (b) *SR*, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d*. Colours represent different node-sizes, horizontal axis is the initial number of particles per node computed from LFN parameters (node-size and particle density) and vertical axis is the temporal average of (1 − *r*) among 19 starting dates. Pearson correlation coefficients *r* were calculated among every pair of node-size/particle-density and its reference pair, which has the same node-size and the highest particle density (5.176 *km*<sup>-2</sup>). Parameters used are described in table 1. The error bars indicate the standard error of the mean and black line represents the function 1/(Initial number of particles per node). The grey dotted line represent the threshold of 0.95 used to identify the recommended number of particles.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Figure 4: Maps of the Local Retention generated with particle density=5.176 km⁻², starting date 15 November 2010, PLD=60 days and 12 m depth. Node size is: (a) 1/2°, (b) 1/4° and (c) 1/8°.
<figure><img src="image_8.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>c)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 5: Maps of the Source Sink strength generated with the same parameters as in Fig. 4 (i.e. particle density=5.176 km⁻², starting date 15 November 2010, PLD=60 days and 12 m depth). Node size is: (a) 1/2°, (b) 1/4° and (c) 1/8°. Values close to zero imply a larval source whereas values close to 1 indicate larval sinks.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>c)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 6: Maps of the Source Sink degree generated with the same parameters as in Fig. 4 (i.e. particle density=5.176 $km^{-2}$, starting date 15 November 2010, PLD=60 days and 12 m depth). Node size is (a) 1/2°, (b) 1/4° and (c) 1/8°. Nodes with values close to zero have a source character, whereas those with values close to 1 indicate a sink character.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 7: Starting-date sensitivity tests for each connectivity metric: (a) LR, (b) SR, (c) $SS_s$ and (d) $SS_d$. Colours represent node-size and PLD, horizontal axis is the time-lag and vertical axis is the temporal average of $r$ among 51 starting dates (1 September to 21 October 2010). Pearson correlation $r$ was calculated among every pair of connectivity metric fields obtained from parameters that only differ in starting dates by the indicated time-lag. The error bars indicate the standard deviation among the set of coefficients obtained.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 8: Spawning-frequency sensitivity tests. Mean relative error (MRE, Eq. (1)) for each connectivity metric: (a) *LR*, (b) *SR*, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d</sub>. Colours represent node-size and PLD as indicated in the legend. The horizontal axis is the period of spawning (inverse of spawning frequency). The vertical axis is the MRE (expressed as a percentage). Lines indicate the MRE of the mean values while error bars represent the MRE of their standard deviations. The reference simulation with respect to which errors are computed is one of daily spawning.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 9: Percentage of discriminated pair of sites (normalized by the number of pairs discriminated using daily spawning frequency) as a function of spawning period for the following metrics: (a) *LR*, (b) *SR*, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d*. Blue curves represent the 0.05 significance level of the Kruskal-Wallis tests; upper and bottom contours correspond to 0.10 and 0.01 significance levels, respectively. Simulations used here are those of the spawning-frequency test (table 1)), with node size of 1/2° and PLD of 60 days. The decreasing trend of the curves indicates that the statistical properties of different sites become more difficult to distinguish for larger spawning period.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 10: Results of PLD sensitivity test for each connectivity metric: (a) LR, (b) SR, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d</sub>. The color represents Pearson correlation among two simulations that only differ by their PLD (in days, given by both axis). Black lines contour the 0.9 level. The temporal averaging concerns 20 starting dates, i.e. 1 September and 1 October over 2001-2010. Other parameters are: particle density of 0.518 km<sup>-2</sup> (i.e. 1600 per node) and node-size 1/2° (see table 1).</figcaption></figure>
# Supplementary Material for “Sensitivity and robustness of larval connectivity diagnostics obtained from Lagrangian Flow Networks”

by Pedro Monroy, Vincent Rossi, Enrico Ser-Giacomi, Cristóbal López and Emilio Hernández-garcía.

## Additional figures

<figure><img src="image_18.png" /><figcaption>Figure SM-1: Number of particles in each node at the end of simulation for node-size=1/8°, starting date = 1 September 2010, initial particle density = 5.176 km⁻² (i.e. 1000 particles at each node), depth = 12 m and PLD = 60 days. Bathymetric contour (black line) at -12 m depth.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

Figure SM-2: Initial particle density sensitivity tests for each connectivity metric: (a) *LR*, (b) *SR*, (c) *SSs* and (d) *SSd*. All nodes in the Mediterranean were used to build this figure (no filtering). Colours represent node-size, X-axis is the initial number of particles per node estimated from LFN parameters (node-size and density particles) and Y-axis is the temporal average of (1−*r*) among 19 starting dates. Pearson correlation *r* was calculated among every pair of node size/density particle and his reference pair, which corresponds to the same node-size and highest initial density particle (5.176 km⁻²). The error bars indicate the standard error of the mean and black line represents the function 1/(Initial number of particles per node).
<figure><img src="image_20.png" /><figcaption>Figure SM-3: Initial particle density sensitivity tests for each connectivity metric: (a) *LR*, (b) *SR*, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d</sub>. This figure was built using only the nodes with land ratio = 0.0 and with longitudes comprised between 1°W and 22°E. Colours indicate node-size, X-axis is the initial number of particles per node estimated from LFN parameters (node-size and density particles) and Y-axis is the temporal average of (1 − *r*) among 19 starting dates. Pearson correlation *r* was calculated among every pair of node size/density particle and his reference pair, which corresponds to the same node-size and highest initial density particle (5.176 km<sup>-2</sup>). The error bars indicate the standard error of the mean and black line represents the function 1/(Initial number of particles per node).</figcaption></figure>
<figure><img src="image_21.png" /></figure>

Figure SM-4: Evolution of connectivity diagnostics when varying the spawning dates over 3 months (1 September 2010 to 30 November 2010) at daily increments (black numbers) for two exemplary nodes: a) node located off Genoa (44°N/8°E) and b) node located over the Ebro shelf (40°N/0°E). Other parameters are kept constant: node size 1/2°, PLD 60 days and density particle 0.518 km⁻². X-axis indicates LR, Y-axis represents SR and the colours represent *SS*<sub>s</sub>. Panel (a) shows relatively stable connectivity proxies when considering successive starting dates (low temporal variability). In contrast panel (b) reveals large and complex fluctuations of connectivity metrics reflecting the variability of ocean currents.

<figure><img src="image_22.png" /></figure>

Figure SM-5: Evolution of connectivity diagnostics when varying the PLD over 5-65 days at daily increments (black numbers) for two exemplary nodes: a) node located over the Ebro shelf (40°N/0°E) and b) node located over the Gulf of Lion shelf (43°N/5°E). Other parameters are kept constant: particle density of 0.518 km⁻² (i.e. 1600 per node) and node-size 1/2°. X-axis indicates LR, y-axis represents SR and the colours represent *SS*<sub>s</sub>. Panel (a) intuitively shows that both retention proxies (LR and SR) decrease when increasing tracking time. Conversely, panel (b) reveals that more complex fluctuations of connectivity proxies can superimpose on top of the overall decreasing trends in retention.
<figure><img src="image_23.png" /></figure>

Figure SM-6: Percentage of discriminated pair of sites (normalized by the number of pairs discriminated with daily spawning frequency) as a function of spawning period for the following connectivity metrics: (a) *LR*, (b) *SR*, (c) *SS*<sub>s</sub> and (d) *SS*<sub>d*. Simulations used here are those corresponding to the spawning frequency test (tab. 1), with node size of 1/2° and PLD of 30 (blue line) and 60 days (red line). Significance level for the Kruskal-Wallis tests is 0.05.