1 Introduction

Spatial patchiness of plankton is a key factor enhancing fishery production and carbon export in the ocean (Okubo and Levin 2001; Franks 2005; Benoit-Bird and McManus 2012; Woodson and Litvin 2015). The formation of plankton patchiness requires some mechanisms for plankton to overcome turbulent diffusion (Abraham 1998). These mechanisms can be broadly classified into two categories: physical and biological. For example, Woodson and Litvin (2015) have shown that convergent currents at ocean fronts create hotspots of plankton, which in turn attract active, larger animals such as fish, sea birds and marine mammals. However, Folt and Burns (1999) argued that physical mechanisms alone cannot fully explain patterns of plankton patchiness. They proposed four biological mechanisms that may be responsible for zooplankton patchiness: diel vertical migration (DVM), predator avoidance, food patches, and mate seeking. Of these, DVM is believed to be “one of the most widespread and powerful biological causes of vertical patchiness” (Folt and Burns 1999).

Zooplankton DVM in the ocean is the greatest animal migration on the Earth in terms of biomass. DVM is believed to be mainly driven by light and can be attributed to a number of mechanisms including avoidance of both predators and UV damage, with the goal of optimizing fitness throughout the life cycle (Cohen and Forward Jr 2016; Thygesen and Patterson 2018). How DVM can induce vertical patchiness is obvious: zooplankton aggregate at certain depths in response to the light stimuli. However, it is not clear whether zooplankton DVM may induce horizontal patchiness.

A pioneering study by Genin et al. (1988) demonstrated that the interactions among zooplankton DVM, seamount, and fish predation can induce horizontal patchiness of zooplankton. This is because the zooplankton above a seamount cannot migrate into deep water, they can be vulnerable to predation by demersal fish. This can produce horizontal “holes” of zooplankton above the seamount. In such a scenario, all the three elements, variable sea bathymetry, predation, and DVM, are essential for the formation of horizontal patchiness. Are all these factors necessary for inducing zooplankton horizontal patchiness? Here, we hypothesize that DVM alone may be able to induce horizontal patchiness.

Our hypothesis is inspired by the insight that the correlation of zooplankton velocities can help to resist diffusion and maintain swarms (Okubo and Anderson 1984). We hypothesize that when zooplankton perform DVM in a water column with vertical gradients of horizontal current velocity, the synchronized vertical positions would lead to synchronized horizontal current velocities that drive zooplankton and subsequently induce horizontal aggregation.

To test this hypothesis, we take advantage of an existing two-dimensional (2D) physical circulation model which can generate strong vertical gradients of current velocity as a result of internal tides. This model was constructed to simulate how the interaction between the Kuroshio and barotropic tidal forcing can generate internal waves at the Izu–Ogasawara Ridge (Masunaga et al. 2019). The Izu–Ogasawara Ridge is a chain of seamounts that stretches over 600 km from the mainland of Japan to the Mariana Trench, separating the Philippine Sea plate to the west from the Pacific plate to the east. This area is strongly affected by the Kuroshio and internal tides generated by the interactions between barotropic tides and seamounts (Masunaga et al. 2017). It is thus an ideal place to study how the interactions between baroclinic flows and zooplankton DVM may induce horizontal patchiness.

In our hypothesis, the vertical gradient of horizontal current velocity is the key to the formation of horizontal patchiness of migrating zooplankton. To test this, we explored two scenarios. The first is that there is no vertical gradient of horizontal velocity of currents (i.e., the water column is unstratified). The second is a realistic stratification case which generates vertical velocity gradients, namely internal waves in this situation. We also test whether the length scale of vertical migration plays a role in affecting horizontal patchiness by creating four zooplankton groups with differing migration amplitude. We also take advantage of the existence of the seamount to examine the effect of uneven bathymetry on the formation of zooplankton horizontal patchiness.

To model plankton patchiness, an individual based modelling approach in the Lagrangian framework has important advantages over the Eulerian framework (Yamazaki and Okubo 1995). Simple biotic interactions on the individual level can induce complex emergent spatial patterns that cannot be easily foreseen (Young et al. 2001). Therefore, we constructed an individual based model of migrating zooplankton and embedded it within the 2D circulation model. We find that both DVM and the seamount can induce horizontal patchiness, via distinct mechanisms.

2 Materials and methods

2.1 2D physical circulation model

To simulate tidal currents and internal waves, we constructed a 2D (horizontal–vertical, or x–z) version of the Stanford Unstructured Nonhydrostatic Terrain-following Adaptive Navier–Stokes Simulator (SUNTANS; Fringer et al. 2006). We compared two cases, one control without stratification and the other with stratification. In the control case, we imposed a constant density with temperature set to 15 °C and salinity set to 30 across the whole model domain, to prevent generation of internal waves. For the stratification case, we computed the full dynamics of currents, temperature and salinity following the Navier–Stokes equation and mass transport equations.

The horizontal model domain spanned from − 500 km (west) to 500 km (east). The horizontal grid spacing was set at 100 m over the ridge and was expanded toward the west and east boundaries to 10 km. The vertical grid spacing was refined near the surface, with a minimum of 3.2 m at the surface increasing to 62.0 m for the deepest grid. The time step was 10 s to maintain numerical stability.

We computed the vertical viscosity and eddy diffusivity using the Mellor–Yamada level 2.5 turbulent closure model (Mellor and Yamada 1982). The horizontal eddy viscosity and diffusivity were set to 0.01 m2 s−1 and 0 m2 s−1, respectively. To generate internal waves, barotropic tidal forcing was imposed at both the western (left) and eastern (right) boundaries with M2 tidal frequency and tidal amplitude of 0.5 m.

The model bathymetry was smoothed from the zonal line along 32.7° N crossing the Izu–Ogasawara ridge, which has the potential to be generalized to any seamount (Fig. 1). To better examine the interaction between DVM and bathymetry in the model, the shallowest summit was lifted to 130 m. To reduce computational costs, the maximum depth was set to 3000 m.

Fig. 1
figure 1

Horizontal velocities after 40 days of simulation with the 2D SUNTANS model. a Unstratified case, b stratified case. White regions indicate the ridge

Full size image

The model was initialized from no motion and was run for 40 days, which was sufficient to generate internal waves. Initial temperature and salinity profiles were obtained from the World Ocean Atlas (WOA) 2009 (Locarnini et al. 2010; Antonov et al. 2010). All the detailed model setup and physical outputs have been reported in Masunaga et al. (2019).

3 Zooplankton diel vertical migration model

3.1 Zooplankton isolume and DVM amplitude

We set four groups of zooplankton with different DVM amplitudes. The first group was a control without DVM. The second to the fourth groups were designed to have different isolumes (preferred light intensity) of 15, 0.04, and 2 × 10–11 W m–2. These values correspond to maximal migrating depths of 50, 149, and 506 m, respectively, assuming a light attenuation coefficient of 0.06 m−1 of seawater. These three migrating zooplankton groups may represent chaetognaths (shallow migrator), copepods (intermediate migrator), and euphausiids (deep migrator) (Genin et al. 1988), respectively. As such, Group 2 would not reach the sea bottom during migration. A minor fraction of Group 3 would encounter the summit of the ridge during migration, whereas a much larger fraction of Group 4 would encounter the ridge during migration. We did not consider other growth and mortality processes of zooplankton, in order to focus solely on the interactions between DVM, currents, and bathymetry.

3.2 Light attenuation with depth

Light was assumed to decrease exponentially with depth following the Beer–Lambert law:

$$ I\left( z \right) = I_{0} e^{{ - k_{{\text{w}}} z}} , $$
(1)

in which I0 was the surface light and kw was the attenuation coefficient (0.06 m−1). I0 was modeled as a sinusoidal function of time of the day between sunrise and sunset (Fig. 2a):

$$ I_{0} = \left\{ {\begin{array}{*{20}l} {I_{{\text{m}}} \sin \left[ {\frac{{\pi \left( {t - t_{r} } \right)}}{{L_{D} }}} \right], \quad {\text{if}}\quad t_{{\text{r}}} \le t \le t_{{\text{r}}} + L_{{\text{D}}} ,} \\ {0, \qquad\qquad\qquad {\text{otherwise}},} \\ \end{array} } \right. $$
(2)

where Im is the maximal light intensity (300 W m–2) at noon on June 1 at 32.7° N latitude. tr is the timing of sunrise (6 a.m.). LD is the day length (14 h).

Fig. 2
figure 2

a Daily variations of light intensity at sea surface. The relationships between the relative zooplankton vertical speed and the distance, b between the current and desired depth and c from surface or bottom boundary

Full size image

3.3 Zooplankton random walk

Motions of zooplankton particles due to advection and diffusion were modeled as a corrected random walk (Visser 1997; Ross and Sharples 2004; Yamazaki et al. 2014):

$$ z_{n + 1} = z_{n} + \left[ {w\left( {z_{n} } \right) + \frac{{\partial K\left( {z_{n} } \right)}}{\partial z}} \right]{\text{d}}t + R\sqrt {\frac{2}{r}K\left( {z_{n} + \frac{{\frac{{\partial K\left( {z_{n} } \right)}}{\partial z}{\text{d}}t}}{2}} \right){\text{d}}t} , $$
(3)

in which Zn+1 and Zn represent the locations along the z axis at n + 1 and n time steps, respectively. \(w\left( {z_{n} } \right) = w_{{\text{b}}} \left( {z_{n} } \right) + w_{{\text{p}}} \left( {z_{n} } \right),\) in which wb represents the swimming speed and wp represents the current velocity, each at depth zn. \(K\left( {z_{n} } \right)\) represents the eddy diffusivity at depth Zn. R is the normalizing factor corresponding to a dimensionless and uniformly distributed random variable with mean zero and variance r, which equals 1/3 for a uniform distribution between − 1 and 1 (Visser 1997; Ross and Sharples 2004). Note that Eq. (3) applies to both horizontal and vertical movements. dt represents the time step (30 s). At each time step, the current velocities and eddy diffusivities were linearly interpolated from the physical outputs, which had a temporal resolution of 30 s.

Vertical migration speed of zooplankton was modeled as a function of the relative distance to the preferred distance of zooplankton (Zp), and the distances to surface (Zs = 0 m) and sea bottom (Zb) (Batchelder et al. 2002). We assumed zooplankton swimming speed would be directly proportional to its distance from its isolume (i.e., the distance between the current position and that of optimal light intensity for a given zooplankter), following Batchelder et al. (2002). To avoid abrupt changes in swimming speed near the boundary (sea surface or bottom), we also assumed that zooplankton swimming speed would be reduced when it was close to sea surface or bottom (Batchelder et al. 2002; Kimmerer et al. 2014). The specific formulations are:

$$ w_{{\text{b}}} = w_{{\text{m}}} \delta \frac{{\left| {Z - Z_{{\text{p}}} } \right|}}{{\left| {Z - Z_{{\text{p}}} } \right| + K_{{\text{p}}} }}\frac{{\left| {Z_{{\text{c}}} } \right|}}{{\left| {Z_{{\text{c}}} } \right| + K_{{\text{c}}} }}, $$
(4a)
$$ \left| {Z_{{\text{c}}} } \right| = \min \left( {\left| {Z - Z_{{\text{s}}} } \right|,\left| {Z - Z_{{\text{b}}} } \right|} \right), $$
(4b)
$$ z_{{\text{p}}} = \frac{{\ln \frac{{I_{0} }}{{I_{i} }}}}{{k_{{\text{w}}} }}, $$
(4c)

where Z is the current vertical location of zooplankton and Ii is the isolume. The maximal vertical swimming speed (wm) was set as 0.035 m s–1 (Batchelder et al. 2002). δ = 1 (upward) if Z is deeper than Zp and δ = − 1 (downward) if otherwise. Kp (= 100 m) and Kc (= 3 m) are the half-saturation constants determining the parabolic function between w and \(\left| {Z - Z_{{\text{p}}} } \right|\) and \(\left| {Z_{{\text{c}}} } \right|\) (Fig. 2).

3.4 Initial and boundary conditions

The model was initialized from midnight (0 a.m.), and a time step of 30 s was used. Two hundred thousand individuals of each group were initialized uniformly within the region from − 200 to 200 km horizontally and from − 200 m (or sea bottom whichever is shallower) to − 5 m vertically. Initial positions were the same for all four groups.

We assumed that zooplankton were passively absorbed onto the boundaries when they were displaced into the boundaries. We also tried the reflective boundary condition in which the zooplankton were reflected back into the water column and found that the different boundary conditions only had slight effects on the model results due to the relatively fast migration speed of zooplankton.

3.5 Analysis of model output

We first examined the patterns of ocean currents and zooplankton distribution after 40 days of simulation to obtain an intuitive view of the distribution patterns. Then we computed the standardized Morisita index (Imst) to quantify the extent of horizontal dispersion (Morisita 1959, 1962; Smith-Gill 1975). To further examine the effects of bathymetry on zooplankton dispersion, we also compared the distributions of bottom depth of the locations of zooplankton individuals after simulation with the initial distributions. Finally, we computed the correlations among horizontal velocities of zooplankton to better understand the relationship between dispersion and correlations of velocity. According to Okubo and Anderson (1984), the dispersion quantified by the variance of particle positions depends on the temporal correlations among particle velocities. If the individual particle velocities are more strongly correlated, they will be less dispersed. In other words, if the particle velocities are not temporally correlated, their positions will gradually diverge with time and vice versa.

To separate the effect of DVM from bathymetry, we separated the horizontal domain into three segments: west (− 200 to − 50 km), middle (− 50 to 50 km), and east (50–200 km). In both the west and east segments, none of the zooplankton individuals would touch the sea bottom in DVM; while in the middle segment, some of the Group 3 and 4 zooplankton would encounter the seamount during vertical migration. Therefore, we assume that in both the west and east segments, the differences of horizontal distributions among different groups would arise purely from the DVM without any bathymetry effect, while in the middle segment, the differences would be caused by both bathymetry and DVM effects.

3.5.1 Standardized Morisita index (I mst)

To calculate Imst, we first calculated the Morisita index (Imor):

$$ I_{{{\text{mor}}}} = N\frac{{\mathop \sum \nolimits_{i = 1}^{N} m_{i}^{2} - \mathop \sum \nolimits_{i = 1}^{N} m_{i} }}{{\left( {\mathop \sum \nolimits_{i = 1}^{N} m_{i} } \right)^{2} - \mathop \sum \nolimits_{i = 1}^{N} m_{i} }}, $$
(5)

in which N is the number of samples. mi is the counts of individuals in the ith sample. The value of Imor ranges from 0 to N and a larger value indicates greater patchiness. Then we normalized Imor between − 1 and 1 to obtain Imst (Smith-Gill 1975):

$$ I_{{{\text{mst}}}} = \left\{ {\begin{array}{*{20}l} {0.5 + \frac{{I_{{{\text{mor}}}} - M_{{{\text{clu}}}} }}{{2\left( {N - M_{{{\text{clu}}}} } \right)}},\quad\; {\text{if}}\quad I_{{{\text{mor}}}} \ge M_{{{\text{clu}}}} > 1} \\ {\frac{{I_{{{\text{mor}}}} - 1}}{{2\left( {M_{{{\text{clu}}}} - 1} \right)}}, \qquad\qquad {\text{if}} \quad M_{{{\text{clu}}}} > I_{{{\text{mor}}}} \ge 1} \\ { - \frac{{I_{{{\text{mor}}}} - 1}}{{2\left( {M_{{{\text{uni}}}} - 1} \right)}},\qquad\quad {\text{ if}}\quad M_{{{\text{uni}}}} < I_{{{\text{mor}}}} < 1} \\ { - 0.5 + \frac{{I_{{{\text{mor}}}} - M_{{{\text{uni}}}} }}{{2M_{{{\text{uni}}}} }}, \quad {\text{ if}}\quad I_{{{\text{mor}}}} < M_{{{\text{uni}}}} < 1 } \\ \end{array} } \right., $$
(6)

in which \(M_{{{\text{clu}}}} = \frac{{\chi_{{\text{L}}}^{2} - N + \mathop \sum \nolimits_{i = 1}^{N} m_{i} }}{{\mathop \sum \nolimits_{i = 1}^{N} m_{i} - 1}}\) and \(M_{{{\text{uni}}}} = \frac{{\chi_{{\text{U}}}^{2} - N + \mathop \sum \nolimits_{i = 1}^{N} m_{i} }}{{\mathop \sum \nolimits_{i = 1}^{N} m_{i} - 1}}\). \(\chi_{{\text{L}}}^{2}\) and \(\chi_{{\text{U}}}^{2}\) are the 2.5% and 97.5% quantiles of \(\chi^{2}\) distribution with degrees of freedom of N − 1. An Imst value between − 0.5 and 0.5 suggests a random distribution insignificantly different from a Poisson distribution. \(I_{{{\text{mst}}}} < - 0.5\) suggests a significantly uniform distribution and \(I_{{{\text{mst}}}} > 0.5\) suggests a clumped distribution significantly different from the random distribution (i.e., patchiness). As such, the usage of \(I_{{{\text{mst}}}}\) instead of \(I_{{{\text{mor}}}}\) has the advantage of not depending on sample size. We used the function ‘dispindmorisita’ in the R package ‘vegan’ to calculate Imst.

We treated each 1 km equidistant horizontal interval as one sample in each segment such that N equaled to 150, 100, and 150 in the west, middle, and east segments, respectively. Then we counted the number of zooplankton particles (mi) within each sample. We also tried other interval lengths as one sample unit and found that the choice of sample unit did not qualitatively affect the results.

3.5.2 Correlations of horizontal velocities

To compute the correlations among horizontal velocities experienced by the four zooplankton groups in the three segments, we randomly sampled two individuals from each group within one 1 km equidistant interval at time zero and calculated the Pearson correlation coefficients of their horizontal velocities during the last 20 days of the simulation.

As our initial assumption that the vertical swimming speeds were all identical within each group was unrealistic, we ran another simulation for zooplankton with variable isolumes in the stratified case, assuming a normal distribution of isolumes, with a standard deviation of 10% of the mean (i.e., a coefficient of variation (CV) of 0.1). This allowed us to evaluate the effect of isolume stochasticity on velocity correlations and horizontal dispersion.

4 Results

4.1 Ocean physics

The two simulation cases generated different patterns of ocean currents. For the unstratified case, horizontal currents were homogeneous vertically and were dominated by the barotropic motion (Fig. 1a). Barotropic tides (sea surface gravity tidal waves) behaved as shallow and long waves without generating the vertical gradient of currents. While the maximum tidal current speeds could reach more than 0.2 m s−1 immediately above the summits, they were much slower elsewhere.

In contrast, for the stratified case, the current distribution showed a strong baroclinic structure (i.e., vertical heterogeneity) generated by internal tides, which has important consequences for particle dispersion as will be shown later. Current speeds were intensified by internal tides and were higher than those of the unstratified control case at most locations (Fig. 1). Although the model tidal forcings added at the boundaries were the same for both configurations (with and without stratification), baroclinic internal tides remained within the study area for a longer time than surface (barotropic) tides. This resulted in stronger currents for the stratified case (the maximal current velocity above the seamounts being 0.53 m s−1 in the stratified case compared to 0.24 m s−1 in the unstratified case). The phase speed of internal tides was roughly 100 times slower than that of surface tides.

4.2 General patterns of zooplankton horizontal dispersion and patchiness

Different ocean currents between the two cases generate different patterns of zooplankton horizontal dispersion. In the unstratified case, the first three zooplankton groups showed similar horizontal variances with initial conditions after 40 days of simulation (Figs. 3a, 4a–d). By contrast, Group 4 showed some aggregations above the trough and on both sides of the ridge (Figs. 3a, 4e). Imst values indicated that although the horizontal distributions of all groups significantly deviated from the random distribution except for Group 1 in the east segment, Imst values were close to the threshold value 0.5 in both east and west segments (Table 1). Imst values in the west segments were slightly higher than in the east segments, which was consistent with the stratified case. This related to the smoother bathymetry slope in the west segments, which generated stronger waves (Fig. 1a). The aggregations were most pronounced in the middle segment where ocean currents were enhanced and some zooplankton encountered the sea bottom, particularly for Group 4.

Fig. 3
figure 3

Two-dimensional distributions of zooplankton individuals after 40 days of simulation. a Unstratified case, b stratified case without isolume stochasticity. The results of the stratified case with isolume stochasticity are qualitatively similar to b and are therefore not shown. Note that the horizontal positions of the particles have been jiggled to facilitate visualization

Full size image
Fig. 4
figure 4

The kernel density distributions and variances (Var, unit: km2) of horizontal locations of four zooplankton groups. a–e unstratified case. f–j stratified case without isolume stochasticity. The results of the stratified case with isolume stochasticity are qualitatively similar to b and are therefore not shown

Full size image
Table 1 Standardized Morisita indices of zooplankton horizontal dispersion
Full size table

In the stratified case, the horizontal distributions were more aggregated than in the unstratified case due to the stronger internal tides (Figs. 3b, 4f–j). Comparing different groups, the passive particles were less clumped than those performing DVM and the zooplankton particles generally showed stronger aggregation with greater DVM amplitudes as indicated by both larger Imst (Table 1) and visual inspections (Fig. 4f–j). Comparing the three segments in the stratified case, Imst were usually larger in the middle segment than in the west and east segments. One exception was that Imst of Group 4 were similar in the west segment compared to the middle segment, which also related to smoother bottom slope that led to the enhanced internal tides in the west segment. In summary, we observed that deep migrators tended to be more clumped and zooplankton were more aggregated when they could encounter sea bottom.

4.3 Bathymetry effect

We observed that deeply migrating zooplankton tended to move away from shallow waters. The kernel density distributions of the bottom depths of the zooplankton positions after 40 days of simulation showed that Group 4 zooplankton clearly avoided shallow waters (< 300 m) in both unstratified and stratified cases compared to initial conditions (Fig. 5d, h). Group 1 zooplankton also tended to move away from shallow waters in the stratified case (Fig. 5e).

Fig. 5
figure 5

The kernel density distributions of the bottom depths of the zooplankton positions before and after simulation for the unstratified case (a–d) and the stratified case without isolume stochasticity (e–h). Deepest depths below 1000 m are not shown to focus on the regions relevant to the DVM range

Full size image

4.4 Correlation of zooplankton velocities

To examine whether horizontal dispersion was associated with correlation of horizontal current velocities, we compared the temporal correlations of horizontal velocities of zooplankton in three segments among the unstratified case and the two stratified cases without and with isolume stochasticity (Fig. 6). The correlations were highly positive (Pearson correlation coefficients near unity) in the unstratified case. The correlations were much more scattered for Group 1 than for other groups in the two stratified cases. However, the correlations were all very high (> 0.99) in other groups, although were slightly lower in the stratified case with isolume stochasticity.

Fig. 6
figure 6

Pearson correlation coefficients of horizontal velocities between two zooplankton individuals randomly selected from 1 km equidistant horizontal intervals in three segments in (a–d) the unstratified case, (e–h) the stratified case without isolume stochasticity (CV = 0), and (i–l) the stratified cases with some isolume stochasticity (CV = 0.1). Note that Y-scales differ

Full size image

5 Discussion

We find that DVM substantially enhances zooplankton horizontal patchiness. Synchronization of horizontal velocities is the key to zooplankton horizontal patchiness, provided that a vertical gradient of horizontal currents exists. This velocity synchronization can be caused by either DVM itself or by the collision of zooplankton with the sea bottom. Below we first discuss the limitations that need to be realized before fully interpreting our results. Then we explain the mechanisms causing velocity synchronization and horizontal patchiness in more detail and further discuss the implications of our findings in a broader context.

5.1 Limitations

The first limitation is the simplified biology in the zooplankton model. In order to focus solely on the effects of DVM on horizontal patchiness, we have neglected any growth, mortality, and other biological behaviors of zooplankton. These factors are potentially important determinants of zooplankton distributions. For example, spatial aggregation may induce greater competition for prey or attract predators (Genin et al. 1988). Particularly, when the deep migrating zooplankton swim out of the euphotic zone, they may experience food limitation (Batchelder et al. 2002). Also, the low temperature at depth may reduce the growth of zooplankton (Loose and Dawidowicz 1994). These effects of low growth and/or high mortality may counterbalance the tendency of migrating zooplankton to aggregate. Currently, we are working to include both bottom-up and top-down factors in the model.

The second limitation is the unique physical setting of the model. This 2D physical model focuses on simulating the internal tides at the Izu–Ogasawara Ridge (Masunaga et al. 2019). It remains to be investigated in detail how the intensity and spatial structure of ocean baroclinic currents induce different patterns of zooplankton horizontal dispersion. It would also be interesting to examine how ocean mixing in a more realistic 3D environment would affect zooplankton dispersion. Although 3D simulations can provide further details of physical processes, the primary internal wave dynamics are similar for both 2D and 3D simulations (see Masunaga et al. 2018, 2019).

The third limitation is the difficulty of completely separating the effects on patchiness of bathymetry and DVM, respectively. One possibility would be to modify the bathymetry to examine the effect of bathymetry on zooplankton patchiness. However, the problem is that the patterns of ocean currents depend on bathymetry. We also tried to compare the patterns of zooplankton patchiness in the three different segments with different bathymetry. However, a similar problem remains: the ocean currents differ among the segments. However, while it is nearly impossible to completely separate these two effects, we can infer from our results that bathymetry and DVM each independently affect zooplankton horizontal patchiness, as discussed below.

5.2 Two mechanisms inducing velocity synchronization and horizontal patchiness

Two mechanisms can be identified contributing to the synchronization of the current velocities experienced by zooplankton and their horizontal patchiness. The first is that DVM synchronizes zooplankton vertical positions and horizontal velocities in a stratified water column. Stratification causes vertical gradients of current velocity (i.e., internal tides). In such a flow field, passive particles at different vertical positions experience different horizontal velocities and therefore tend to disperse horizontally over time. By contrast, the vertical positions of migrating zooplankton tend to track their isolumes and are synchronized by the light cue. As such, the horizontal velocities they experience tend to be correlated (i.e., synchronized), which induces horizontal patchiness (Okubo and Anderson 1984). Adding stochasticity of migrating amplitudes decreases correlation of horizontal velocities (Fig. 6) and also reduces the extent of horizontal patchiness (Table 1), further confirming the linkage between velocity synchronization and horizontal dispersion.

Zooplankton horizontal patchiness increases with migrating amplitude (Table 1; Fig. 3). This trend is reflected by the differences of correlation coefficients of horizontal velocities between Group 1 and three migrating groups (Fig. 6). However, the differences of correlation coefficients are too small to explain the differences of Imst observed in the three groups of migrators. While part of the enhanced synchronization of deep migrators might be accounted for by interactions with sea bottom as discussed below, it can also arise from other factors such as increased migration speed and depth-dependent horizontal velocity patterns. This is evidenced by comparing the Morisita indices in the east and west segments where no zooplankton would encounter the sea bottom but horizontal patchiness still increases with migration amplitude (Table 1). It appears that the deep migrators tend to be trapped within some localities in the west segment due to the convergence nature of the internal tides they experienced.

The second mechanism is the interaction between bathymetry and DVM. Vertically migrating zooplankton tend to pursue their optimal environment (light in this case), which is a deterministic force that drives them toward the deepest environment they can reach when light level increases. If they do not encounter the sea bottom, their horizontal locations will be subject to only horizontal advection and diffusion. However, if the migrating zooplankton encounter the sea bottom, their horizontal positions will be affected if the sea bottom is not flat. During the period when surface light level is increasing, migrating zooplankton will consistently move towards deeper waters along the oblique sea bottom to pursue their optimal light level, like a ball rolling downhill. In contrast, when the surface light decreases and the zooplankton ascend, they will leave the bottom and their horizontal positions will not be affected by bathymetry. As a consequence, zooplankton will not be able to return to their original horizontal position if they encounter an uneven sea bottom. The overall net effect will be a tendency for the migrating zooplankton to move toward locations with deeper bottom depth (Fig. 3b), so long as it is shallower than the depth of their optimal light level. This point is also valid for Group 1 because although they could not active migrate, they were also subject to vertical diffusion, which also made them move away from shallow waters.

It is worth noting that the two above mechanisms are relatively independent of each other. The uneven sea bottom alone can induce zooplankton horizontal aggregation without the existence of baroclinic currents (internal tides/waves). In the unstratified case with vertically homogeneous horizontal current velocities, Group 4 zooplankton still form some aggregations atop the ridge due to the avoidance of shallow waters (Table 1; Fig. 3a). DVM alone can also cause horizontal aggregation even in the deep open ocean as long as baroclinic currents exist. This is supported by the aggregations of zooplankton in the west and east segments where zooplankton do not encounter the sea bottom (Table 1).

5.3 Implications for secondary production and fisheries

Spatial patchiness is an important factor contributing to productivity (Woodson and Litvin 2015). Since water column stratification and seamounts can promote zooplankton patchiness, these two factors should also enhance secondary production and fisheries. It has been observed that aggregations of migrating zooplankton trapped within the trough of seamounts may attract predatory fish and become feeding hotspots for marine fish near seamounts (Genin et al.1988, 1994). Seamounts have been already reported as important sites for enhanced productivity and fisheries because of enhanced nutrient supply via N2 fixation or upwelling (Hasegawa et al. 2008, 2009; Shiozaki et al. 2014). This study provides evidence that seamounts can also enhance secondary production by promoting zooplankton patchiness.

6 Conclusion

We find that zooplankton DVM can induce horizontal patchiness in a water column with vertically heterogeneous horizontal ocean currents (i.e., baroclinic currents). Larger DVM amplitudes lead to less dispersion and greater patchiness. The main underlying mechanisms are that: (1) DVM tends to synchronize the horizontal current velocities experienced by zooplankton, and (2) deep migrators tend to move away from shallow waters when they encounter the sea bottom. To our knowledge, these seemingly simple mechanisms have not been mentioned in the literature.

Our findings suggest that DVM may be a ubiquitous cause of zooplankton horizontal patchiness in the ocean, which enhances fishery production and carbon export (Okubo and Levin 2001; Franks 2005; Benoit-Bird and McManus 2012; Woodson and Litvin 2015). Our modelling results suggest that the isolume is a key trait of zooplankton and should be frequently measured in field studies. Our results also generate a working hypothesis for future field programs: the patchiness of zooplankton species should be positively associated with their DVM amplitudes and isolumes.