Uncovering marine connectivity through sea surface temperature

The δ-MAPS analysis is performed onto monthly mean SST anomalies from the Mediterranean Sea Physical Reanalysis (CMEMS MED-Physics²⁵) over the period 1987–2017. The advantage of using a reanalysis resides in the availability of a velocity field consistent with the SSTs that allows us to confirm the coupling between network domains and ocean currents within the euphotic layer.

Validation of the δ-MAPS framework

The proposed ecoregionalization is first applied to the 2007–2010 period, when the domains, representing ecoregions, can be compared to those identified by Ref.² using Lagrangian methods. The details of this validation are reported below and relevant figures can be found in the Supplementary Information.

The 2007–2010 ecoregions in Figure S1 are consistent with the ones derived in Ref.² through computationally intensive simulations. The name of each domain corresponds with those used in Ref.² to ease comparison. It is worthwhile remarking that this work and the one of Ref.² not only use very different methods to define connectivity, but also different data sources. Our study uses velocity and SST output fields from CMEMS MED Physics reanalysis, while Ref.² uses the configuration PSY2V3 of the operational system MERCATOR OCEAN with a resolution of 8 km in the horizontal downscaled to a connectivity grid of 50 × 50 km. The data assimilation and clustering algorithms are different and Ref.² employs a cut-off in addition to the clustering grid downscaling. These differences unavoidably translate into slightly different shapes and patterns of the domains inferred. For example, the D + V area in panel (a) of Figure S1 is effectively two separate ecoregions in Ref.², in which the Messina Strait is not resolved at the connectivity grid level. However, this separation appears inconsistent with the surface kinetic energy (K.E.) of panel (b) in Figure S1, computed from the horizontal currents, e.g. zonal (u) and meridional (v) velocity components, as K.E. = 1/2 |V|² where |V|= (u² + v²)^0.5. Indeed, there is no clear separation between the regions north and south at Messina Strait in our dataset. Having detailed this example and acknowledged that some differences should be expected, the overall basin eco-regionalization using δ-MAPS is consistent with that in Ref.². The spatial accuracy is enough to well separate the main ecological areas, despite small-scale differences (i.e. some km, due to resolution choices).

By and large, the SST anomaly domains in Figure S1 are bounded by ocean currents, in agreement with Ref.². This is due to the dominance of advective forcing by ocean currents on the SSTs at equatorial and mid latitudes, on monthly timescales and spatial scales of few hundreds kilometers²⁴. This link, which is foundational to the proposed methodology, is further quantified as follows: First, we calculate the surface K.E. per unit mass averaged over the time slot of interest; second, we select the points in the validation period (2007–2010) that exceed the 50^th percentile of surface K.E. computed for the entire basin over 1987–2017 (e.g. 0.004 m²/s²); third, we compute the domain-boundary matrix augmented by 1 grid point in each direction; finally, we count which fraction of the domain boundaries computed in the boundary matrix overlaps with the K.E. fronts (above the 50th percentile threshold). The fraction obtained is high and equal to 0.73, and remains elevated when increasing the threshold to the 60th percentile (0.66). This procedure was repeated for all the time slots with Δ = 7 years used next in this study, obtaining high and very stable values in each case (mean ± variance = 0.73 ± 0.01 for the 50th percentile threshold, and 0.65 ± 0.01 for the 60th percentile threshold).

Additionally, the correlation between the surface K.E. and K.E. at 50 m, 100 m, and 150 m over the whole 1987–2017 period (Figure S2) remains positive and significant, with coefficients for the whole domain (field mean c.c.  ± variance) of 0.83 ± 0.04 at 50 m, 0.68 ± 0.05 at 100 m, and 0.54 ± 0.06 at 150 m, indicating that the link extends to the whole euphotic layer.

Mediterranean Sea ecoregions: long-term changes

The space-averaged (e.g. averaged on the whole basin) SSTs over the 1987–2017 period are characterized by a linear warming trend of about 0.04 °C per year, stronger in the eastern portion of the basin (Figure S3 in Supplementary Information). Over the same period, the K.E. per unit mass is characterized by different trends over decadal or quasi-decadal periods (Fig. 2, shown for surface only but the trend extends similarly to 50 m and 100 m depths) and no clear east–west contrast. A positive trend is found in the first part of the curve (1987–2001, 2.3 × 10^–4 m²/s² per year, green line in figure), followed by a central decade without statistically significant changes (2001–2010, blue line), and a steep negative trend afterward (2010–2017, – 4.1 10^–4 m²/s² per year, red line). We refer to 1987–2001, 2001–2010 and 2010–2017, as the UP, MAX and DOWN periods. The dynamical changes associated with the strengthening and weakening of ocean currents are hypothesized to coincide with a reshaping of the sub-basin ecoregions and reciprocal connectivity. The ecoregionalization inference is therefore performed considering time slots of varying length, so that yr_end = yr_ini + Δ with yr_ini = y₀ + n, n = 0,1,…,N, where y₀ is the initial year of the dataset (1987) and N is the total number of time slots, each of duration Δ years, between 6 and 8. Time slots overlapping by more than one year among different trends periods are excluded. The choice of Δ = 7 years represents the best trade-off for having enough time slots to quantify the evolution of ecoregions and a sufficiently large number of data points in each time slot for statistical inference. We will focus on this case, but results are verified also for the other Δ values (see Supplementary Information).

Strength maps for three representative time slots are presented in Fig. 3a,c,e while maps of domain strengths for all Δ = 7 time slots can be found in Figure S4. The mean surface kinetic energy averaged within each timeslot is next compared to the number of ecoregions in corresponding timeslots. The fragmentation level, or the total number of ecoregions, and the mean surface kinetic energy content are highly correlated (Figure S5b in Supplementary Information), with a Pearson’s coefficient of 0.79 for the whole Mediterranean Sea, and 0.8 (0.65) for the eastern (western) basin. The fact that time slots are not independent does not invalidate the analysis, and a large correlation (c.c = 0.73) is retained even when using four non-overlapping time slots. A higher fragmentation occurs whenever the upper ocean layer is more energetic, and this relationship is robust to changes of Δ (see Supplementary Information). The domain strength is next compared to the mean K.E. content. For each timeslot, the domain strength is spatially averaged over the eastern and western basin separately. The correlations between the averaged strengths and the corresponding time slot mean surface K.E. values, both varying as the time slots change, are then calculated for eastern and western basins separately. No linkage is found in the western basin, but a strong anticorrelation describes the relationship in the eastern Mediterranean (c.c. − 0.74). This anticorrelation remains high (− 0.73) also when the eastern basin strengths are related to the whole basin surface K.E. averaged over each timeslot.

We hypothesize that the amount of K.E. associated with semi-permanent jets, currents or large mesoscale eddies, grouped here together and named KE fronts, can be used as an indicator of their role as connectivity modulators. We identify KE fronts applying a pattern recognition algorithm on the K.E. fields for each time slot. The resulting pictures are processed by an image segmentation technique, based on K-means clustering, to separate the K.E. in four clusters of increasing energy content. The maximum-intensity group is selected as indicators for KE fronts and the number of pixels contained in each cluster is counted and used to estimate the size or abundance of each one. The maximum-intensity cluster well represents the energy-containing structures as measured by the correlation between the mean surface K.E. content in each time slot and the pixels within the corresponding cluster (c.c. > 0.99). The more pixels reside within each cluster, the larger the KE fronts-populated areas that this cluster approximates. This estimation is carried out for the whole basin, and separately in the eastern and western parts. The number of pixels is then correlated to the number of inferred ecoregions for the whole Mediterranean (c.c. = 0.81), and for eastern (c.c. = 0.81) and western (c.c. = 0.69) basins. Figure S6 in the Supplementary Information compares the clustering maps of a low energy time slot (1987–1993, in panel (a)) and a higher one (2004–2010 in panel (b)), for the whole Mediterranean Sea for the maximum cluster. The number of ecoregions is highly correlated with the KE fronts everywhere and especially in the eastern Mediterranean Sea. The higher level of fragmentation found in the MAX period is thus associated with more abundant and/or larger surface KE fronts, acting as eco-dynamical barriers.

To further strengthen this assessment, we consider that energy fronts can act as modulators for SSTa-derived domains. The ecoregionalization over a certain time slot characterizes that time range in one single ecoregion-map but stems from data known at several time points (i.e. monthly SSTa in our case). The resulting domains account therefore for the inherent physical variability of the system over time. A higher (lower) ecoregions fragmentation may therefore by associated with dynamical fronts occurring at different times and not necessarily in the same place, over a certain time range. If this is plausible, we expect to count more (less) occurrences of higher energy in broad areas where the domains are more (less) fragmented. For each time slot, the number of occurrences of a front in each pixel is therefore counted. Specifically, having defined a front as a K.E. realization above the 50th percentile of the overall (1987–2017) time varying surface K.E., we count how many times a front appears in the considered time slot at each pixel. In Figure S7 pixels are colored according to the number of occurrences in each time slot. The result is consistent with the domain fragmentation evolution. The higher fragmentation occurring in timeslots from 2001 to 2010 in the eastern basin is associated with more frequent fronts. Similar considerations hold for the other sub-basins, clearly distinguishing low energy periods from higher ones.

Mediterranean ecoregions connectivity networks

Changes in functional networks or connectivity among ecoregions can be assessed by comparing a network from each energy period (UP: 1987–1993, MAX: 2004–2010 and DOWN: 2011–2017) (Fig. 3b,d,f for the eastern basin and Figure S8 in the Supplementary Information for the western basin).

In 1987–1993 the western basin was characterized by a high mean positive correlation of 0.73, with a strong, non-directional connectivity among the Tyrrhenian and Ligurian-Algero Provençal domains. In 2004–2010 the connectivity was overall weaker, and in particular reduced among Tyrrhenian waters. The connectivity between the Balearic domain (Bal) and the Tyrrhenian ones was also reduced. In 2011–2017 the connectivity was mostly recovered, especially in Tyrrhenian waters. In this period, the Algero-Provençal domain separated from the Ligurian Sea (Lig), enforcing its connectivity with the Balearic and the Alboran ecoregions.

In the eastern basin we focus our attention on the ecoregion immediately offshore the Suez Canal (Fig. 3), the major anthropogenic corridor for the introduction of non-indigenous marine species in the Mediterranean Sea, the so-called Lessepsian immigrants³². According to δ-MAPS, connectivity from the domain surrounding Suez was high in the first decade, decreased approaching MAX, remained small until about 2010–2011 with fewer statistically significant links, and increased again in the more recent time slot considered. During the UP and DOWN periods, the strongest connections were with the eastern Levantine (domain N), followed by that with the Aegean, Ionian and Tunisian Seas. During UP the connectivity extended to the Provençal and Algerian Seas, in the western basin, while in DOWN these links were absent and replaced by a connection with the Adriatic Sea.

The 1987–1993 and 2011–2017 periods, while not too dissimilar in energy levels, differed indeed for the phase of the Ionian-Adriatic Bimodal Oscillating System or BiOS^33,34. The BiOS is a mode of variability characterized by a decadal reversal of the Northern Ionian Gyre (NIG) from cyclonic to anticyclonic, and vice versa. In its anticyclonic spinning the NIG deviates the inflowing Modified Atlantic Water (MAW) from the Sicily Channel towards the northern Ionian, entering the Adriatic Sea and decreasing its salinity and temperature. This prevents a portion of the MAW from reaching the Levantine basin, and enhances the outflow of Levantine waters into the western basin, along a pathway that follows the African coastline. The anticyclonic NIG co-occurs with higher concentrations of Atlantic and Western Mediterranean organisms in the Adriatic Sea. When the NIG is cyclonic, on the other hand, Levantine waters enter the Adriatic Sea, whereas the MAW preferably flows toward the Levantine³⁵ and Lessepsian migrations influence the Adriatic Sea at various latitudes, affecting also phytoplankton phenology^33,36,37. The corresponding regions and connectivity networks in the two opposite NIG periods are detailed in Figure S9 in the Supplementary Information.

Source: Ecology - nature.com