Chlorophyll shading reduces zooplankton diel migration depth in a high-resolution physical–biogeochemical model

Zooplankton diel vertical migration (DVM) is critical to ocean ecosystem dynamics and biogeochemical cycles, by supplying food and injecting carbon into the mesopelagic ocean (200–800 m). The deeper the zooplankton migrate, the longer the carbon is sequestered away from the atmosphere and the deeper the ecosystems they feed. Sparse observations show variations in migration depths over a wide range of temporal and spatial scales. A major challenge, however, is to understand the biological and physical mechanisms controlling this variability, which is critical for assessing impacts on ecosystem and carbon dynamics. Here, we introduce a migrating zooplankton model for medium and large zooplankton that explicitly resolves diel migration trajectories and biogeochemical fluxes. This model is integrated into the MOM6-COBALTv2 ocean physical–biogeochemical model and is applied in an idealized high-resolution (9.4 km) configuration of the North Atlantic. The model skillfully reproduces observed North Atlantic migrating zooplankton biomass and DVM patterns. Evaluation of the mechanisms controlling zooplankton migration depth reveals that chlorophyll shading decreases zooplankton migration depths by 60 m in the subpolar gyre compared with the subtropical gyre, with pronounced seasonal variations linked to the spring bloom. Fine-scale spatial effects (<100 km) linked to eddy and frontal dynamics can either offset or reinforce the large-scale effect by up to 100 m. This could imply that, for phytoplankton-rich regions and filaments, which represent a major source of exportable carbon for migrating zooplankton, a high chlorophyll content contributes to reducing zooplankton migration depth and carbon sequestration time.

1 Introduction

Diel vertical migration (DVM) of zooplankton is the largest synchronized movement of organisms on Earth, with major implications for ocean ecosystem dynamics and biogeochemical cycles (Hays et al., 1997; Steinberg et al., 2002; Bianchi et al., 2013b). Every day at dawn, after spending the night feeding, half of the world's zooplankton biomass leaves the surface ocean for deeper waters, where it remains before swimming back to surface waters at dusk (Klevjer et al., 2016). Migrating zooplankton shape biogeochemical cycles by consuming organic matter in the surface waters and respiring or excreting some of it deeper down, sequestering carbon away from the atmosphere and potentially accounting for 10 %–30 % of the biological carbon pump (Steinberg et al., 2000; Boyd et al., 2019; Aumont et al., 2018; Archibald et al., 2019; Nowicki et al., 2022). This efficient vertical transport of organic matter by DVM influences the depth at which nutrients are regenerated (Steinberg et al., 2002; Longhurst and Glen Harrison, 1988; Bianchi et al., 2013b, 2014) and intensifies respiration and oxygen depletion at the upper margin of oxygen minimum zones (around 200–500 m deep; Bianchi et al., 2013a; Aumont et al., 2018). DVM also contributes to structuring trophic interactions by providing energy to mesopelagic ecosystems (Kelly et al., 2019), triggering cascade migration of lower trophic levels (Bollens et al., 2011) and altering the foraging performance of visual predators such as fish and marine mammals by hiding in the dark (Benoit-Bird and Moline, 2021; Chambault et al., 2024).

Zooplankton migration depth is a key factor in understanding the impact of DVM on carbon and oxygen cycles as well as ecosystem dynamics. The deeper migrating zooplankton egest, excrete or respire material from the surface, the deeper nutrients are regenerated and oxygen is consumed and the longer the carbon they release is sequestered away from the atmosphere (Boyd et al., 2019). The fecal pellets they egest can also serve as a food source for ecosystems at that depth and below (Kelly et al., 2019). Migration depth observations vary by several hundred meters over multiple spatial and temporal scales. For example, at fine spatial scales, Powell and Ohman (2015) showed that zooplankton migrate 200 m deeper on one side of a 20 km front than on the other side. At the basin scale, Bianchi et al. (2013a) showed that zooplankton swam about 200 m deeper in the subtropical gyre than in the North Atlantic subpolar gyre. Temporally, Omand et al. (2021) showed that migration depth can fluctuate by up to 50 m within a day, while Hobbs et al. (2021) documented the existence of a seasonal cycle in zooplankton migration depth. However, the mechanisms controlling zooplankton diel vertical migration patterns and their biogeochemical consequences are still poorly quantified (Bandara et al., 2021).

Light is often considered the main cue that determines zooplankton migration timing and depth (Brierley, 2014). The most common hypothesis, known as the “preferendum hypothesis”, states that migrating organisms modify their depth when irradiance changes so as to remain at a preferred light level or isolume (Ewald, 1910; Michael, 1911; Russell, 1927). This behavior keeps them hidden from visual predators (Zaret and Suffern, 1976). Numerous studies supporting this hypothesis have shown that factors modifying irradiance at the ocean surface, such as sunlight, moonlight, cloud cover, ice cover or eclipses (Strömberg et al., 2002; van Haren and Compton, 2013; Last et al., 2016; Omand et al., 2021; Flores et al., 2023), or in the water column, such as water transparency and chlorophyll content, could modify the migration depth (Dickson, 1972; Powell and Ohman, 2015; Hobbs et al., 2021). Migration timing can be modified further by other exogenous factors such as predation pressure (Ohman et al., 1983; Lass and Spaak, 2003; Bollens and Frost, 1989) or food availability (Pijanowska and Dawidowicz, 1987; Beklioglu et al., 2008) or by endogenous ones such as circadian rhythms (Harris, 1963). Similarly, migration depth can also be modified by exogenous parameters such as low oxygen concentrations (Bianchi et al., 2013a), temperature or salinity (Kimmerer et al., 1998; Glaholt et al., 2016), together with endogenous ones such as zooplankton size (Ohman and Romagnan, 2016). However, the limited number of observations makes it challenging to disentangle and quantify the individual contributions of these different factors to migration depth spatial and temporal variations.

Although numerical models can provide further insights, modeling diurnal vertical migration poses several challenges. First, migration is affected by several physical, chemical and biological variables (e.g., light, oxygen concentration, food availability and predation exposure), requiring models that capture these complex processes as well as a realistic representation of zooplankton physiology that is key to biogeochemical impacts (e.g., egestion, excretion, respiration, growth and mortality rates). Second, diel vertical migration is a rapid process where zooplankton reach migration velocities ranging from 3 to 15 cm s⁻¹ (Bianchi and Mislan, 2016), and it is modulated at various temporal and spatial scales. Modeling such fast zooplankton vertical migration over a wide period and ocean area at high spatial and temporal resolution is therefore limited by the computational cost. Existing models fail to capture this complexity. One-dimensional vertical models explicitly represent zooplankton physiology and migration trajectory based on isolumes (Bianchi et al., 2013b; Nocera et al., 2020) or take into account exposure to food and predators to maximize fitness (Pinti et al., 2019), but they do not represent ocean circulation. Conversely, three-dimensional global biogeochemical models capturing ocean circulation represent vertical migration in parameterized form, where the biogeochemical fluxes created are distributed around a migration depth determined by an isolume or oxygen threshold (Aumont et al., 2018; Archibald et al., 2019; Nowicki et al., 2022). In this case, migration trajectories are not explicitly modeled, and the seasonal cycle is not always represented. Furthermore, none of these modeling approaches captures the effect of fine-scale dynamics (< 100 km), such as fronts and eddies, on zooplankton migration patterns. To improve our understanding of the impact of DVM on biogeochemical and ecosystem dynamics, we need to develop models explicitly representing the migration trajectories and physiology of migrating zooplankton, together with their interactions with the environment over a wide range of space scales and timescales.

In this study, we develop a migrating zooplankton module and fully integrate it into a coupled physical–biogeochemical model of the North Atlantic. This model includes a range of biological and physical mechanisms that affect zooplankton migration patterns and physiology over multiple spatial and temporal scales. We use an idealized double-gyre high-resolution configuration of the North Atlantic, built using an ocean circulation model (Adcroft et al., 2019) coupled with a biogeochemical model (Stock et al., 2020). Our idealized configuration replicates the characteristic biophysical dynamics of the subtropical and subpolar gyres, their seasonal cycle and their fine-scale dynamics. Our migrating zooplankton module includes two migrating zooplankton size classes (medium and large) and explicitly represents their migration trajectories, which are controlled by light, oxygen concentration and size and prey distribution. Zooplankton physiology is also explicitly modeled and depends on environmental factors (e.g., temperature) and individual factors (e.g., zooplankton feeding activity or swimming). We show that our zooplankton migration model reproduces the contrasts in the biomass, timing and migration depth of migrating zooplankton observed in the North Atlantic. We identify and quantify the mechanisms responsible for modulating the migration depth of zooplankton in the North Atlantic (surface irradiance, chlorophyll shading and ocean transport). We explore their variation across seasons, biomes (subtropical and subpolar) and spatial scales (gyre scale and fine scales) and show that chlorophyll shading is the dominant process structuring migration depth at all these scales.

2 Materials and methods

We develop an idealized double-gyre ocean model reproducing the North Atlantic biophysical ocean dynamics (Sect. 2.1). Our model is coupled to the biogeochemical module (COBALTv2) and extended in this study to integrate two vertically migrating zooplankton (COBALTv2-DVM, Sect. 2.3). Our new module represents realistic migration trajectories (Sect. 2.3.1), visual predation (Sect. 2.3.2) and zooplankton physiology (Sect. 2.3.3). We also develop a framework to quantify the processes modulating zooplankton migration depth (e.g., irradiance seasonality, the chlorophyll shading effect at large and fine scales, or vertical transport; see Sect. 2.4).

2.1 Idealized double-gyre experiment

2.1.1 Configuration

We use an idealized double gyre reproducing biophysical dynamics characteristic of the North Atlantic Ocean, including a low-productivity subtropical gyre and a high-productivity subpolar gyre separated by a jet analogous to the Gulf Stream, as well as a deep winter convection region (Fig. 1). Idealized double-gyre models are useful tools for examining the influence of gyre circulation and eddies on ventilation (Marshall et al., 2002; Radko and Marshall, 2003; Henning and Vallis, 2004) as well as phytoplankton and biological carbon pump dynamics (Lévy et al., 2015; Resplandy et al., 2012, 2019; Couespel et al., 2021).

Here, the double-gyre model was built using the Modular Ocean Model version 6 (MOM6, Adcroft et al., 2019). The domain is a square of 3570 km side length, 4000 m depth, 9 km horizontal resolution, 2–5 m vertical resolution between 0 and 100 m depth, 10–100 m between 100 and 1000 m depth and up to 250 m below. The configuration is centered around 40° N and under the β-plane approximation ($f_{\mathrm{0}}=\mathrm{9.35}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹, $\mathit{\beta }=\mathrm{1.75}\times {\mathrm{10}}^{\mathrm{11}}$ m⁻¹ s⁻¹). The four vertical boundaries are solid impermeable walls preventing any lateral exchanges. This physical model is coupled to the Carbon, Ocean Biogeochemistry and Lower Trophics biogeochemical module version 2 (COBALTv2, Stock et al., 2020). COBALTv2 includes 33 biogeochemical tracers, such as the main nutrients (e.g., nitrogen, phosphorus, iron, or silica), and a planktonic ecosystem. The original version of COBALTv2 included three phytoplankton and three non-migrating zooplankton, to which we added two migrating zooplankton (see Sect. 2.3). Vertical movement was simulated with an implicit solution of the advection–diffusion equation generalized to handle upward and downward swimming. The resulting tridiagonal matrix was solved using the algorithm described in Press (2007). Each year, the double gyre is forced by the same idealized zonally uniform atmospheric climatology detailed below. Further details on the model parameters are in Appendix A (Table A1).

Figure 1 Double-gyre setup: 10 June snapshot of (a) sea surface height (top) and depth sections of temperature at 20° N (front) and nitrate concentration at −15° E (right), (b) sea surface temperature, (c) mixed-layer depth, (d) eddy kinetic energy and (e) surface chlorophyll concentration.

2.1.2 Forcings

The ocean surface is forced with an idealized wind field (Fig. A1a in Appendix A). The zonal wind stress follows a sinusoidal profile that is minimal at the northern and southern boundaries of the domain (−0.1 N m⁻²), maximal at 40° N (0.1 N m⁻²) and varying by 10 % in amplitude according to a seasonal cycle (e.g., higher values in December). The meridional wind stress is zero. The radiative forcing is derived from the ERA5 atmospheric reanalysis (Hersbach et al., 2020) averaged over the 1979–2022 period (Fig. A1g). The temporal resolution of ERA5 (1 h) accurately resolves the daily cycle of irradiance required for vertical migration modeling. The forcing is calculated as the zonal mean over a thin band of 2° longitude located in the middle of the North Atlantic (longitude −35 to −33° E, latitude 20 to 60° N), which prevents the daily light cycle from being skewed by the time difference between the eastern and western Atlantic basin. This radiative forcing, like all the other surface forcings, is temporally synchronous across the whole basin on daily timescales and linearly interpolated every 5 min between the available forcing values. This resolution enables our model to accurately represent the mean daily light cycle, although it may not perfectly resolve the crepuscular period, when zooplankton begin to migrate. The other physical atmospheric forcings (pressure, temperature, precipitation and humidity) are derived from the JRA55-do atmospheric reanalysis (Tsujino et al., 2018) with a 3 h time resolution averaged over the period 1958–2020 (Fig. A1b–f). They are calculated as the zonal mean over a wider region covering most of the North Atlantic Ocean (longitude −60 to −20° E, latitude 20 to 60° N). Precipitation from JRA55-do is adjusted to ensure that the same mass of precipitation is received by the double gyre as by the North Atlantic Ocean, and a mass restoration flux is applied at the surface to ensure its conservation (Fig. A1e). Iron and lithogenic material deposition (Fig. A1h–i) is calculated as the zonal mean of the last 30 years of a GFDL-ESM4 preindustrial control simulation spanning the time from 1850 to 2014 (Dunne et al., 2020) over the same band as the physical forcings.

2.1.3 Initial state and spinup

Temperature and salinity as well as nitrate, oxygen, phosphate and silica concentrations are initialized uniformly and zonally using the World Ocean Atlas 2018 zonal mean (Boyer et al., 2018a) over the region delimited by −60 to −20° E and 20 to 60° N. Similarly, alkalinity and dissolved inorganic carbon are initialized from GLODAPv2 (Olsen et al., 2020), and the rest of the biogeochemical tracers are initialized from the GFDL-ESM4 simulation described above, all using zonally averaged fields in the same region. The model is first spun up for 1000 years at a coarse resolution (85 km) and then run for 25 years at an eddy-permitting resolution (9 km, Table A1). This paper focuses on the last 5 years at 9 km resolution.

2.2 Observations

To evaluate the biophysical dynamics simulated in the model, we use multiple observational data products. The nutrient concentration, temperature and salinity are taken from the World Ocean Atlas 2018 (Garcia et al., 2019, Fig. A2), the mixed-layer depth from a compilation of hydrographic sections between 1941 and 2002 in de Boyer Montégut et al. (2004), the surface chlorophyll concentration from a combination of satellite observations between 1997 and 2020 in Sathyendranath et al. (2019) and the net primary productivity from Moderate Resolution Imaging Spectroradiometer (MODIS) measurements between 2002 and 2023 processed by three algorithms (Eppley, CbPM and CAFE; Behrenfeld and Falkowski, 1997; Westberry et al., 2008; Silsbe et al., 2016).

The biomass of migrating zooplankton is evaluated using estimates derived from lidar observations between 2007 and 2019 (Behrenfeld et al., 2019) and net samples at the Bermuda Atlantic Time Series (BATS) station between 1994 and 2010 (Steinberg et al., 2012). We also use estimates of migration depth and time spent at depth by migrating zooplankton derived from acoustic Doppler current profiler (ADCP) measurements taken between 1990 and 2010 and compiled by Bianchi and Mislan (2016). Note that these observations are particularly rare and challenging to collect and may be subject to measurement method bias. Acoustic methods can be biased towards highly sound-scattering animals (Flagg and Smith, 1989), net measurements do not capture animals smaller than the mesh or larger ones able to avoid the net and lidar data can be biased towards larger plankton and miss biomass capture of smaller plankton.

2.3 Migrating zooplankton model

We developed the COBALTv2-DVM module by adding two migrating zooplankton to the pre-existing COBALTv2 biogeochemical module (Stock et al., 2020). These migrating zooplankton are the same as the medium (0.2–2 mm) and large (2–20 mm) zooplankton of COBALTv2 in terms of prey preferences and density-dependent higher-trophic-level predator losses. However, they perform diel vertical migration (Sect. 2.3.1), enabling them to escape visual predation (Sect. 2.3.2). They also have a compartmentalized physiology that temporally decouples ingestion, digestion, respiration and growth (Sect. 2.3.3) to realistically represent the matter fluxes which differ when zooplankton feed at the surface and when they dive, as highlighted by the work of Bianchi et al. (2013b).

2.3.1 Migration model

The migration trajectory of zooplankton is controlled by the light intensity, oxygen concentration and vertical distribution of their prey (Fig. 2a). Migration is assumed to be solely vertical and is controlled by the vertical migration velocity (w).

Figure 2 (a) Top of the panel – mean daily surface irradiance cycle (black line). Bottom of the panel – mean daily cycle of migrating zooplankton biomass (medium: green color fill, large: red color fill), migrating zooplankton isolume depth (medium: green line, isolume = 10⁻³ W m⁻², large: red line, isolume = 10⁻⁴ W m⁻²) and irradiance (blue color fill). (b) Time spent at depth by migrating zooplankton as a function of day length in observations (cross) and models (distribution of the model points shown by grey color fill).

During the day, migrating zooplankton dive following an isolume, i.e., a preferred light intensity level (I_iso; medium I_iso = 10⁻³ W m⁻², large $I_{\mathrm{iso}}={\mathrm{10}}^{-\mathrm{4}}$ W m⁻²). As a result, large zooplankton dive deeper than medium zooplankton. The migration speed w(z), positive downward, is close to the maximum swimming speed (w₀) away from the isolume and starts to decrease, according to a Monod equation, when it approaches about 50 m from its isolume:

$$\begin{array}{}\text{(1)}& w\left(z\right)={\mathit{\eta }}_{\mathrm{iso}}\left(z\right)w_{\mathrm{0}}{\displaystyle \frac{\left|\frac{\mathrm{1}}{K_{\mathrm{Blue}}}\ln \right(\frac{I_{\mathrm{iso}}}{I\left(z\right)}\left)\right|}{K_{\mathrm{d}z}+\left|\frac{\mathrm{1}}{K_{\mathrm{Blue}}}\ln \right(\frac{I_{\mathrm{iso}}}{I\left(z\right)}\left)\right|}},\end{array}$$

with K_blue=0.0232 m⁻¹ the light attenuation coefficient for optically pure seawater (Manizza et al., 2005), K_dz=50 m the half-saturation distance, I(z) the irradiance at depth z and w₀ the maximum migration velocity determined from the observations (medium w₀=6 cm s⁻¹ and large w₀=8 cm s⁻¹; Fig. A7). η_iso(z) is 1 or −1, depending on the zooplankton position (z) relative to the isolume depth and their prey distribution, and controls the migration direction. When zooplankton are above the isolume depth (z<z_iso), η_iso(z)=1 and zooplankton swim downwards. When the zooplankton are below the isolume, e.g., because the isolume rises at dusk or is absent at night, migrating zooplankton return to the surface to feed. To reproduce foraging behavior, migrating zooplankton distributions are redistributed according to the distribution of their prey and η_iso(z) can be 1 or −1, depending on

$$\begin{array}{}\text{(2)}& \mathrm{ZooExcess}\left(z\right)=\underset{\mathrm{0}}{\overset{z}{\int }}\left({\displaystyle \frac{Z\left(z^{\prime }\right)}{\stackrel{\mathrm{‾}}{Z}}}-{\displaystyle \frac{P\left(z^{\prime }\right)}{\stackrel{\mathrm{‾}}{P}}}\right)\mathrm{d}{z}^{\prime },\end{array}$$

where Z(z^′) and P(z^′) are the concentrations of zooplankton and their prey at depth z^′. $\stackrel{\mathrm{‾}}{Z}$ and $\stackrel{\mathrm{‾}}{P}$ are their mean concentrations over the entire water column. At a given depth z₀, if ZooExcess(z)>0, the zooplankton are in relative excess with respect to their prey between the surface and this depth. Therefore, zooplankton dive (η_iso(z)=1) to restore the balance. If ZooExcess(z)<0, the zooplankton are in deficit relative to their prey. Therefore, zooplankton rise to fill this imbalance (${\mathit{\eta }}_{\mathrm{iso}}\left(z\right)=-\mathrm{1}$). Note that, if the zooplankton below the isolume rise and reach it, the zooplankton remain at the isolume depth, since their migration velocity becomes zero ($w(z=z_{\mathrm{iso}})=\mathrm{0}$).

The model also includes a migration depth limitation by an oxygen threshold set to 10 µmol O₂ kg⁻¹, similarly to Aumont et al. (2018). However, this limitation is not utilized in our simulation as oxygen concentrations are higher than 10 µmol O₂ kg⁻¹ and therefore do not limit migration. The sensitivity of the results to the migration model parameters (w₀, I_iso and K_dz) is assessed by varying their values by ±50 % (see Fig. B3 in Appendix B).

2.3.2 Visual predation model

In the ocean, zooplankton migrate vertically to escape visual hunting predators. To represent this competitive advantage of migrating zooplankton over non-migrating zooplankton, we implemented the visual predation model derived by Bianchi et al. (2013b) in COBALTv2-DVM. In this model, a fraction of the predation rate of top predators (I_HP) on zooplankton (Z), corresponding to visual predation, increases with irradiance (I(z), Fig. A8):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{I}_{\mathrm{HP}}& =I_{\mathrm{HP}}^{\max }\cdot {\displaystyle \frac{Z}{K_Z+Z}}\\ & \cdot \left((\mathrm{1}-\mathit{\alpha })+\mathit{\alpha }{\displaystyle \frac{I\left(z\right)}{\left(\frac{K_Z}{K_Z+Z}\right)\cdot K_{\mathrm{irr}}+I\left(z\right)}}\right),\end{array}\end{array}$$

where $I_{\mathrm{HP}}^{max}$ is the maximum predation rate of top predators that varies with temperature and oxygen concentration, K_Z=1.25 µmolN kg⁻¹ is the half-saturation constant for prey density dependence, α=0.9 is the maximum predation fraction due to visual predators and K_irr = 10⁻¹ W m⁻² is the half-saturation response for light limitation. At the surface during the day, irradiance is sufficient (I(z)≫K_irr) and does not limit visual predation: $I_{\mathrm{HP}}\left(z\right)\approx I_{\mathrm{HP}}^{\max }\frac{Z}{K_Z+Z}$. During the night or at depth, the absence of light (I(z)≪K_irr) inhibits visual predation and limits the predation rate to 10 % of its maximum value: $I_{\mathrm{HP}}\left(z\right)\approx (\mathrm{1}-\mathit{\alpha })I_{\mathrm{HP}}^{\max }\frac{Z}{K_Z+Z}$.

2.3.3 Zooplankton physiological model

To reproduce realistic biogeochemical fluxes produced by zooplankton digestion, respiration and death, we developed a six-compartment model. The model, inspired by that of Bianchi et al. (2013b), was extended to represent phosphate, iron and silica dynamics as well as the thermal dependence of the gut-clearing rate. The model represents nitrogen (N^gut), phosphorus (P^gut), iron (Fe^gut) and silica (Si^gut) in the zooplankton gut, together with nitrogen metabolites (N^metab) and nitrogen body biomass (N^body, Fig. 3).

Figure 3 Migrating zooplankton model COBALTv2-DVM: (a) simplified and (b) detailed schemes of biogeochemical fluxes and stocks. The gut content (blue box) is filled by ingestion and cleared by egestion, excretion and assimilation, occurring between 10 min and 3 h after ingestion. The metabolite pool (light-red box) is filled through assimilation and used for zooplankton growth and respiration, occurring throughout the day. The biomass (green box) accumulates through growth and declines through mortality.

First, zooplankton fill their guts with the nitrogen, phosphorus, iron and silica contained in their prey. The gut content (X^gut=N^gut, P^gut, Fe^gut or Si^gut) is cleared by egestion and assimilation following an e-folding timescale ranging from 15 min to 3 h, depending on water temperature (Fig. A6):

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial X^{\mathrm{gut}}}{\partial t}}={\mathrm{ingestion}}_X-k^{\mathrm{clear}}\cdot X^{\mathrm{gut}},\text{(5)}& {\displaystyle }{k}^{\mathrm{clear}}=k_{\mathrm{0}}^{\mathrm{clear}}+k_T^{\mathrm{clear}}\cdot (T-T_{\mathrm{0}}),\end{array}$$

where T₀=0 °C and the coefficients $k_{\mathrm{0}}^{\mathrm{clear}}=\mathrm{8}$ d⁻¹ and $k_T^{\mathrm{clear}}=\mathrm{4.32}$ °C⁻¹ d⁻¹ are derived from a linear regression on a set of digestion rate observations for different temperatures (Irigoien, 1998). The ingestion rate (ingestion_X) depends on prey and zooplankton abundance and increases exponentially with temperature (see Stock et al., 2014).

Following the same ratios as the non-migrating zooplankton, 30 % of the nitrogen and phosphorus consumed is egested as organic matter (medium: 70 % of the egestion as particulate sinking detritus and 30 % in dissolved form; large: 100 % of the egestion is detritus). The remaining 70 % of the consumption is assimilated (assimilation_N|P) as metabolites at a ratio of N : P = 1 : 16, allowing the metabolic content to be modeled with a single tracer (N^metab). Nitrogen and phosphorus in excess compared to this ratio are egested in dissolved inorganic form (NH₄ and PO₄, respectively). Silica and iron are egested totally, 30 % as organic matter (medium: 70 % as particulate sinking detritus and 30 % in dissolved form; large: 100 % as detritus), and the rest is egested as dissolved inorganic matter. These ratios are chosen to be the same as those used by Stock et al. (2014, 2020) for non-migrating zooplankton.

Metabolites are consumed by metabolic reactions (metabolism), i.e., the sum of anabolic (growth through organic biomass synthesis) and catabolic (respiration and production of energy necessary for the functioning of the organism) reactions according to an e-folding timescale of 24 h (k^resp = 1 d⁻¹):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial N^{\mathrm{metab}}}{\partial t}}={\mathrm{assimilation}}_{\mathrm{N}|\mathrm{P}}-\mathrm{metabolism},\text{(7)}& {\displaystyle }\mathrm{metabolism}=\mathrm{growth}+\mathrm{respiration}=k^{\mathrm{resp}}\cdot N^{\mathrm{metab}}.\end{array}$$

The respiration flux is the sum of three terms. The first term is the respiration of metabolites to maintain the vital functions (respiration_basal), the second term covers the energetic needs of assimilation and digestion (respiration_feeding) and the last term is linked to swimming (respiration_swim):

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}\mathrm{respiration}& ={\mathrm{respiration}}_{\mathrm{basal}}+{\mathrm{respiration}}_{\mathrm{feeding}}\\ & +{\mathrm{respiration}}_{\mathrm{swim}},\end{array}\text{(9)}& {\displaystyle }{\mathrm{respiration}}_{\mathrm{basal}}={\mathit{\mu }}_{\mathrm{basal}}\cdot N^{\mathrm{body}},\text{(10)}& {\displaystyle }{\mathrm{respiration}}_{\mathrm{feeding}}={\mathit{\gamma }}^{\mathrm{resp}}\cdot k^{\mathrm{resp}}\cdot N^{\mathrm{metab}},\text{(11)}& {\displaystyle }{\mathrm{respiration}}_{\mathrm{swim}}={\displaystyle \frac{w}{w^{\mathrm{ref}}}}\cdot {\mathit{\mu }}_{\mathrm{basal}}\cdot N^{\mathrm{body}},\end{array}$$

where μ_basal is the basal respiration rate of zooplankton that increases exponentially with temperature (see Stock et al., 2014), γ^resp is the fraction of metabolic flux used to cover assimilation and digestion processes, w is the migration velocity of zooplankton and w^ref is the reference velocity for adjusting the energy cost of swimming. The linear relationship between respiration (respiration_swim) and swimming velocity was proposed by Torres and Childress (1983) and previously used in the model developed by Bianchi et al. (2013b). Biomass production by anabolism (growth) is equal to the difference between the metabolic (metabolism) and catabolic (respiration) fluxes. If this difference is positive, the biomass increases. If it is negative, some of the zooplankton do not have the metabolic resources to cover their catabolic needs. This part of the zooplankton dies, and its biomass is transformed into sinking detritus.

Most of the model parameters are derived directly from the literature (such as the digestion rate $k_{\mathrm{0}}^{\mathrm{clear}}$ and its temperature sensitivity $k_T^{\mathrm{clear}}$) or were already part of the pre-existing COBALTv2 model. Four parameters (basal respiration rate μ_basal, grazing rate, reference swimming speed w^ref and maximum predation fraction due to visual predators α) were adjusted using 14 sensitivity experiments (see Fig. A5). We selected the parameter values to match the observed fraction and biomass of migrating zooplankton (see Figs. 5h and A5). Note that this complex physiological parameterization is only applied to migrating zooplankton in order to separate fluxes temporally and therefore spatially. Without migration, this parameterization would give similar results to non-migrating zooplankton, on a daily average.

2.4 Migration depth analysis framework

In this section, we present a framework to investigate which processes control zooplankton migration depth. First, we isolate the effects of surface irradiance, chlorophyll shading and ocean transport (Sect. 2.4.1). Then, we separate their contributions according to different spatial and temporal scales (Sect. 2.4.2).

2.4.1 Migration depth decomposition

We decompose the migration depth (z_mig) as the sum of three terms (Fig. 4) controlled by irradiance in transparent water (z_iso,tw), variations in chlorophyll shading (Δz_chl) and variations in ocean circulation vertical transport (Δz_circ):

$$\begin{array}{}\text{(12)}& z_{\mathrm{mig}}=z_{\mathrm{iso},\mathrm{tw}}+\mathrm{\Delta }{z}_{\mathrm{chl}}+\mathrm{\Delta }{z}_{\mathrm{circ}}.\end{array}$$

z_iso,tw is the isolume depth in transparent water. It represents the depth at which the isolume would be if only the effects of light absorption by water molecules were taken into account. In our model, the coefficients of light attenuation by water are constants (k_red=0.225 m⁻¹ and $k_{\mathrm{blue}/\mathrm{green}}$ = 0.0232 m⁻¹) (Manizza et al., 2005), so z_iso,tw only depends on the ocean surface irradiance.

Δz_chl is the chlorophyll shading effect. It is calculated as the difference between the modeled depth of the isolume and its depth in transparent water ($\mathrm{\Delta }{z}_{\mathrm{chl}}=z_{\mathrm{iso}}-z_{\mathrm{iso},\mathrm{tw}}$). Δz_chl quantifies the reduction in isolume depth due to light absorption by chlorophyll.

Δz_circ is the ocean circulation effect. It is calculated as the difference between the migration depth of zooplankton and their isolume depth at zenith, i.e., when the isolume is deepest ($\mathrm{\Delta }{z}_{\mathrm{circ}}=z_{\mathrm{mig}}-z_{\mathrm{iso}}$). Δz_circ quantifies variations in migration depth caused by vertical physical transport (advection and mixing), which can move zooplankton away from their isolume faster than zooplankton can respond by swimming.

Figure 4 Framework used to identify the processes controlling zooplankton migration depth. The migration depth (z_mig, red line) is the sum of the isolume depth in transparent water (z_iso,tw, yellow line), the chlorophyll shading effect (Δz_chl, green arrow) and the vertical ocean transport effect (Δz_circ, pink arrow). z_iso is the isolume depth targeted by migrating zooplankton. See the Methods section (Sect. 2.4) for details.

2.4.2 Separation of spatial scales

The chlorophyll content and vertical transport, which control migration depth, respond to large-scale dynamics (> 100 km) but are also modulated by fine-scale dynamics (< 100 km) such as mesoscale eddies or submesoscale fronts. To isolate these two scales of variability, we perform a spatial Reynolds decomposition. We spatially filter Δz_chl and Δz_circ using a Gaussian kernel with a standard deviation of 100 km and write them as the sum of a large-scale average (〈Δz_chl〉,〈Δz_circ〉) and a fine-scale anomaly ($z_{\mathrm{chl}}^{\prime },z_{\mathrm{circ}}^{\prime }$) relative to this average:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{\Delta }{z}_{\mathrm{chl}}=\langle \mathrm{\Delta }{z}_{\mathrm{chl}}\rangle +\mathrm{\Delta }{z}_{\mathrm{chl}}^{\prime },\text{(14)}& {\displaystyle }\mathrm{\Delta }{z}_{\mathrm{circ}}=\langle \mathrm{\Delta }{z}_{\mathrm{circ}}\rangle +\mathrm{\Delta }{z}_{\mathrm{circ}}^{\prime }.\end{array}$$

In contrast, the isolume depth in transparent water (z_iso,tw), which depends on surface irradiance, does not vary on fine spatial scales in our simulation because the irradiance forcing is zonally uniform. Consequently, the fine-scale term is almost zero:

$$\begin{array}{}\text{(15)}& z_{\mathrm{iso},\mathrm{tw}}=\langle z_{\mathrm{iso},\mathrm{tw}}\rangle +z_{\mathrm{iso},\mathrm{tw}}^{\prime }\approx \langle z_{\mathrm{iso},\mathrm{tw}}\rangle .\end{array}$$

To summarize, this framework provides a way of disentangling the complexity of the mechanisms modulating the migration depth:

$$\begin{array}{}\text{(16)}& z_{\mathrm{mig}}=z_{\mathrm{iso},\mathrm{tw}}+\langle \mathrm{\Delta }{z}_{\mathrm{chl}}\rangle +\mathrm{\Delta }{z}_{\mathrm{chl}}^{\prime }+\langle \mathrm{\Delta }{z}_{\mathrm{circ}}\rangle +\mathrm{\Delta }{z}_{\mathrm{circ}}^{\prime }.\end{array}$$

We isolate the irradiance contribution and distinguish the spatial contributions of chlorophyll shading and vertical transport at large and fine scales. This framework does not include the limitation of migration depth by oxygen, as its concentration is high and does not limit vertical migration in our simulation.

3 Results

3.1 North Atlantic biophysical dynamics

3.1.1 Contrasted biogeochemical biomes

The double-gyre model reproduces the characteristic physical and biogeochemical dynamics of four biomes of the North Atlantic: a low-productivity subtropical gyre, a high-productivity subpolar gyre, an energetic jet in between with intense mesoscale eddy activity and a deep winter convection region in the northern part of the domain (Figs. 1, 5 and A2). The subtropical, subpolar and jet biomes are defined by thresholds of the annual surface chlorophyll concentration (0.15 and 0.35 mg m⁻³), as used in previous studies (Resplandy et al., 2012, 2019). The convecting region is identified as a subregion of the subpolar region, using the 200 m threshold on the annual mean mixed-layer depth. Our main conclusions are not sensitive to these thresholds. Note that the wind creates an upwelling along the wall at the southern boundary and stimulates ocean productivity, with annual mean chlorophyll concentrations higher than 0.35 mg m⁻³ (Fig. 5d). This upwelling, representative of the equatorial dynamics occurring further south in nature, is excluded from the analysis using the chlorophyll thresholds given above.

Figure 5 Map of the annual mean (a, b) mixed-layer depth (MLD), (d, e) surface chlorophyll concentration and (h, i) migrating zooplankton biomass in (a, d, h) the double-gyre model and (b, e, i) the observation-based data products. The subtropical, jet and subpolar biomes are delimited by the solid black lines (mean surface chlorophyll concentration = 0.15 and 0.35 mg m⁻³), while the subpolar and convective biomes are delimited by a dashed line (mean MLD = 200 m; see the Methods section). Mean seasonal cycle of (c) MLD, (f) net primary productivity and (j) migrating zooplankton biomass in the subtropical biome (orange) and the subpolar biome (blue) in the double-gyre model (solid lines) and observation-based products (dashed lines). The shading and error bars correspond to 2 standard deviations around the average. The orange star in panel (i) indicates the Bermuda Atlantic Time Series (BATS) study. See Sect. 2.2 for the observational product details.

The mean biophysical conditions and seasonal cycle simulated within the biomes are similar to those observed in the ocean. In the observations and the model, the subtropical anticyclonic gyre (elevated sea surface height, Fig. 1a) is oligotrophic at the surface, with nitrogen concentrations averaging 0.29 and 0.28 µmol N kg⁻¹, respectively, and phosphorus concentrations averaging 0.03 and 0.06 µmol P kg⁻¹ in the first 100 m (Fig. A2a–d). The model also reproduces the observed stratification throughout the year, with a mixed layer varying between 30 and 60 m seasonally (Fig. 5a–c), strongly limiting the supply of nutrients to the surface. Consequently, net primary productivity in the subtropical gyre is low (between 170 and 250 mg C m⁻² d⁻¹ on average in the model and observation-based estimates) and varies little seasonally (Fig. 5d–f). In contrast, the subpolar cyclonic gyre (depressed sea surface height) is a nutrient-rich environment. Nitrate concentrations range from 8 to 17 µmol N kg⁻¹ in the top 100 m and phosphate concentrations from 0.3 to 0.8 µmol P kg⁻¹ (Figs. 1a and A2a–d) in the observations and the model. This is partly due to high seasonality, with a modeled mixed layer averaging 140 m in winter (125 m in the observations, Fig. 5c), bringing nutrients to the surface. Consequently, net primary productivity in the subpolar gyre is high and seasonally marked (Fig. 5f). It increases from January onwards, peaks in May (920 mg C m⁻² d⁻¹ in the model, 650–930 mg C m⁻² d⁻¹ in the observations) and falls back below 500 mg C m⁻² d⁻¹ at the end of summer.

Finally, at the northern boundary, the wind creates a downwelling and a deeply convecting water column characterized by winter mixed layers exceeding 1000 m (Fig. 5a). These deep convection events can be considered representative of ocean dynamics occurring further north in the Atlantic Ocean, such as in the Irminger and Labrador seas (Fig. 5b). These convection events bring a high concentration of nutrients to the surface (> 12 µmol N kg⁻¹ and > 0.7 µmol P kg⁻¹; Figs. 1a and A2a–d) and sustain a biome productivity of 420 mg C m⁻² d⁻¹ on the annual average, reaching up to 1250 mg C m⁻² d⁻¹ in spring.

These large-scale dynamics are locally modulated by fine-scale structures such as eddies and fronts (Fig. 1b–e). The jet hosts the most energetic fine-scale dynamics (Figs. 1d, A3), because it is unstable and therefore induces eddy formation. Many eddies are present in the subtropical gyre, since the horizontal resolution of the model is sufficient for resolving mesoscale dynamics (> 100 km). Conversely, in the subpolar gyre, the resolution is not sufficient for fully resolving these eddies, leading to less coherent and more filamentary turbulent structures. As we will show in this study, these features introduce variability into the gyre and seasonal dynamics and are key to understanding zooplankton migration depth variability.

3.1.2 North Atlantic migrating zooplankton dynamics

The model reproduces the regional and latitudinal contrasts in the biomass and fraction of migrating zooplankton observed in the North Atlantic (Figs. 5h–j and A5). In the subpolar biome, the migrating zooplankton biomass and seasonal variations are high in both the model and lidar estimates (1±0.5 g m⁻², Fig. 5c), reflecting the high productivity and seasonality of this biome (Fig. 5d–f). Migrating biomass is minimal in late winter (0.69±0.34 g m⁻²), accumulates in spring and peaks in May (1.17±0.57 g m⁻²). Biomass remains high throughout summer and fall (1.11±0.54 g m⁻²) before declining in November. In contrast, the subtropical biome is characterized by lower biomass and low seasonality in both the model and observations (0.18±0.12 g m⁻², Fig. 5h–j), reflecting the low productivity and seasonality of the region (Fig. 5d–f). Migrating zooplankton account for about half of the medium- and large-sized zooplankton in the model in both biomes, which is in line with observations from the COPEPOD database (Aumont et al., 2018, Fig. A5).

Figure 6 (a) Latitudinal distribution of migration depth. The blue dots are observations derived from acoustic Doppler current profiler (ADCP) observations and the blue line their average per 5° latitude bins. The red line is the zonal and temporal average of the model. The shading corresponds to 2 standard deviations around the average. (b) Migration depth derived from ADCP observations in the North Atlantic used in panel (a).

The model also reproduces the timing of zooplankton migration. The ADCP observations compiled by Bianchi and Mislan (2016) show a linear relationship between the length of the day and the time spent at depth by migrating zooplankton (Fig. 2b). The longer the days, the longer migrating zooplankton remain at depth. This relationship is reproduced in the model. Migrating zooplankton spend between 9 and 18 h at depth, depending on the length of the day. Since the day length is less variable over the year at low latitudes than at high latitudes, the same applies to the time spent at depth by migrating zooplankton. In the subpolar biome, migrating zooplankton stay for averages at depth of 10 h in winter and 18 h in summer, while in the subtropical biome they stay for averages at depth of 11 h in winter and 15 h in summer.

Finally, the model reproduces the contrasts in migration depths observed in the North Atlantic. Zooplankton migration depths vary both across and within biomes. ADCP observations show that zooplankton dive about 100 m deeper in the subtropical biome, where the average migration depth is 480 m, than in the subpolar biome, where it is confined to 370 m (Fig. 6a, b). The model reproduces these large-scale contrasts with average migration depths of 480 m in the subtropical biome and 380 m in the subpolar biome (Fig. 6a). The jet region, around 40° N, shows intermediate migration depths of around 340 m in both the observations and our model. The observations also show migration depth variability within biomes up to 100 to 200 m around the biome average, which is captured relatively well by the model (shading in Fig. 6a). The model simulates migration depth and variability in the convective biome, which are higher than those observed at 55–60° N. The simulated mean migration depth is over 450 m, whereas the observations show depths of around 420 m, and the variability in the model is 3 times greater than in the observations. However, these deep convection events occur further north in nature, as in the Irminger Sea, where migrations have been observed below 500 m (Fig. 6b). Furthermore, in these deep convection regions, some zooplankton replace their daily migratory behavior with seasonal migratory behavior for part of the year (see the discussion in Sect. 4).

3.2 Controls of zooplankton migration depth

The observed migration depth and its variability result from the combination of processes varying on different space scales and timescales, which cannot be disentangled from the limited number of observations available. In the following, we therefore use the model, which reproduces the observed patterns, to investigate the mechanisms controlling the high spatial contrast between the subtropical and subpolar gyres, the variability introduced at fine scales by eddies and filaments (> 100 km) and the influence of seasonality. We identify chlorophyll shading as a major mechanism in regulating migration depth, although the contribution of ocean transport in the convective region and the seasonality of surface irradiance are significant. Only results for medium migrating zooplankton are presented. Results for large migrating zooplankton are similar, but where differences arise, they are indicated as such and presented in the Appendix (Figs. B1 and B2).

3.2.1 Mechanisms modulating zooplankton migration depth

Zooplankton migration depth is modulated by three mechanisms: surface irradiance, chlorophyll shading and ocean vertical transport (see the Methods section, Sect. 2.4). We illustrate the influence of these three effects using model snapshots and a depth–latitude section across the different biomes early in the spring bloom (1 March, Fig. 7). First, migrating zooplankton target a specific isolume whose depth depends on surface irradiance intensity. As surface irradiance increases, light penetrates further into the water column, resulting in a deeper isolume. We track the isolume depth as if the water were transparent (z_iso,tw) to isolate this irradiance effect. For example, Fig. 7 shows that, on 1 March, z_iso,tw is slightly deeper in the subtropical gyre than in the subpolar gyre, because surface irradiance at zenith is stronger at low latitude than at high latitude (yellow line and arrow in Fig. 7d).

Figure 7 Snapshot on 1 March of the (a) biomass-weighted medium zooplankton migration depth at zenith, (b) water column chlorophyll content and (c) mixed-layer depth. The blue line corresponds to the location of the cross section shown in panel (d). (d) Vertical cross section of the chlorophyll (green) and medium zooplankton (red) concentrations on 1 March, at 30° E, between 27 and 57° N. The depth–latitude section also shows the transparent water isolume depth (yellow line), isolume depth (green line) and migration depth (red line) at zenith. The transparent water isolume depth is a baseline controlled by irradiance intensity (yellow arrows), the difference between the isolume and transparent water isolume depths is caused by chlorophyll shading (green arrows) and the difference between migration depth and isolume depth is the effect of ocean vertical transport (red arrows).

Seawater contains chlorophyll, which modifies light penetration in the water column by changing its opacity. As chlorophyll content increases, water opacity and chlorophyll shading increase, resulting in a shallower isolume (z_iso, Fig. 7d, green line). Shallow migration regions are characterized by a high chlorophyll content and deep migration regions by a low chlorophyll content. This chlorophyll shading effect (Δz_chl, green arrow) acts on large scales between the subtropical and subpolar biomes but also at the scale of eddies and filaments, where migration depths can vary by up to 200 m between points that are only 20 km apart, with a chlorophyll content that differs by more than 60 mg m⁻² (Fig. 7a, b).

Finally, migrating zooplankton can be transported away from their isolume when ocean vertical transport associated with advection and/or turbulent mixing is high. For example, in the convective biome, vertical ocean transport (Δz_circ, red arrow) carries migrating zooplankton more than 200 m below their isolume (Fig. 7d). This occurs during deep convection events where the mixed layer exceeds 1000 m (Fig. 7c).

3.2.2 Controls of migration depth spatial variability at large scales

We generalize our analysis by evaluating the average annual effect of irradiance, chlorophyll shading and vertical transport on migration depth across the four biomes (Fig. 8). We find that the spatial contrast in migration depth between the different biomes is mainly due to differences in chlorophyll shading. On an annual average, the effect of irradiance (z_iso,tw) is similar across the four biomes (yellow bars in Fig. 8). Without other factors, irradiance in transparent waters would lead to migration down to 480 m in the subtropical biomes where irradiance is highest and 440 m in the convective biome where irradiance is lowest.

Figure 8 Mechanisms controlling the migration depth of medium zooplankton on an annual average in the (a) subtropical, (b) jet, (c) subpolar and (d) convective biomes. The decomposition includes the effects of large-scale variations in irradiance (z_iso,tw, yellow bar), chlorophyll shading (〈Δz_chl〉, green bar) and vertical ocean transport (〈Δz_circ〉, red bar), together with their seasonal range of variability (black hatch). The effect of transport is insignificant (light-red bar) when it is smaller than the vertical resolution of the model grid at that depth.

Large-scale chlorophyll shading (〈Δz_chl〉) reinforces the weak latitudinal contrast set by irradiance, reducing migration depth by about 100 m on annual average in the less productive subtropical gyre, while it reduces it by approximately 140–165 m in the more productive jet region and subpolar gyre and by up to 180 m in the convective region (Fig. 8a–c, dark-green bars; see Fig. 5d–f for the annual mean chlorophyll). In the convective biome, the effect of ocean vertical transport (〈Δz_circ〉) partially offsets that of chlorophyll shading. Migrating zooplankton are transported 100 m below their isolume on an annual average due to deep convection events. In other regions, the effect of transport is not significant (〈Δz_circ〉 lower than the model vertical grid resolution at the depth of migration; Fig. 8a–c). The migration depth of large zooplankton is 100 m deeper than that of medium zooplankton due to irradiance, since they follow a darker isolume. However, the modulation of their migration depth by chlorophyll shading and transport is similar to that of medium zooplankton (Fig. B1). We note that the chlorophyll shading effect is not sensitive to the choice of parameters in the migration model (e.g., target isolume or migration velocity; Fig. B3).

3.2.3 Controls of migration depth seasonal variability

Migration depth is modulated on the seasonal timescale in response to variations in irradiance, bloom dynamics (i.e., chlorophyll absorption) and winter convection (Fig. 9). Seasonal irradiance variations deepen the migration depth in summer (deep z_iso,tw) and decrease it in winter (shallow z_iso,tw) compared to its average annual effect (Fig. 9, yellow line). The amplitude of seasonal variations in the irradiance effect increases with latitude (subtropical: 25 m, jet: 40 m, subpolar: 57 m and convective: 85 m) due to greater variations in the solar incidence angle.

Large-scale chlorophyll shading considerably reduces migration depth during the spring phytoplankton bloom when the chlorophyll content is highest (Fig. 9, green line). In productive regions, this reduction in the migration depth is most pronounced between February and May and intensifies with latitude due to a more pronounced and prolonged phytoplankton bloom. The seasonal reduction in migration depth attributable to chlorophyll shading reaches 190 m in the jet region, 235 m in the subpolar gyre and 300 m in the convective region (Fig. 9b–d). In contrast, there are no seasonal variations in the chlorophyll shading in the subtropical gyre (Fig. 9a). Large-scale vertical transport deepens the migration depth in the convective region during convective events between January and April, down to −340 m at its peak (Fig. 9, red line). This effect is negligible at other times and in other regions.

Figure 9 Mechanisms controlling the medium zooplankton migration depth seasonal (solid lines) and fine-scale (< 100 km, shading envelope) variability in the (a) subtropical biome, (b) jet, (c) subpolar and (d) convective biomes for the irradiance (z_iso,tw in yellow), chlorophyll shading (〈Δz_chl〉 and $\mathrm{\Delta }{z}_{\mathrm{chl}}^{\prime }$ in green) and vertical transport (〈Δz_circ〉 and $\mathrm{\Delta }{z}_{\mathrm{circ}}^{\prime }$ in red) effects. Contributions are relative to the baseline depth set by mean annual irradiance (isolume in transparent water, decreasing from 520 m in the subpolar gyre to 480 m in the convective region; see Fig. 8).

3.2.4 Controls on migration depth spatial variability at fine scales

Within biomes, patchiness in phytoplankton and mixed-layer depths introduces strong variations in the chlorophyll shading and transport effects ($\mathrm{\Delta }{z}_{\mathrm{chl}}^{\prime }$ and $\mathrm{\Delta }{z}_{\mathrm{circ}}^{\prime }$), and it modulates the migration depth at fine spatial scales, as shown by the snapshots from 1 March (Fig. 7a–c). The amplitude of this fine-scale variability depends on the biome and time of year, with stronger fine-scale variability found in late winter and spring in the jet, subpolar and convective biomes (Fig. 9). On average, the chlorophyll fine-scale effects have a larger impact on migration depth than the fine-scale effects of transport. The only exception is in the convective region, where chlorophyll and transport fine-scale effects are similar.

During the spring bloom, fine-scale variations in chlorophyll and the associated shading (Fig. 9, green filling) modulate the migration depth by ±67 m in the jet, ±75 m in the subpolar gyre and ±100 m in the convective region (Fig. 9b–d, green shading) and are smaller during the rest of the year (less than ±30 m). In other words, this fine-scale effect can create migration depth contrasts of up to 200 m between two spots 10–20 km apart during the spring bloom in productive regions, reinforcing or offsetting the large-scale chlorophyll shading effect. This fine-scale chlorophyll shading effect is however much weaker (±20 m) and shows little seasonal variations in the subtropical gyre waters containing less chlorophyll (Fig. 9a).

In the convective region, the fine-scale effect of transport (Fig. 9, red filling) is high from January to April and modulates migration depths up to ±135 m, due to local changes in the mixed-layer depth (Fig. 7c). During the rest of the year and in other regions, the fine-scale ocean transport effect is not significant. Thus, although these fine-scale effects are small on an annual average (less than ±40 m), they can significantly modulate the migration depth from late winter to spring.

4 Conclusions and discussions

In this study, we have developed a mechanistic zooplankton vertical migration model, which we integrated into a large-scale ocean biogeochemical model and implemented in an idealized, high-resolution double-gyre physical framework. This model reproduces the large-scale contrasts in zooplankton migration depth observed between the North Atlantic subtropical and subpolar biomes (480 m vs. 380 m) (Bianchi et al., 2013b). It also captures seasonal variations similar to those observed by Hobbs et al. (2021), with migration depth decreasing by around 100 m during the phytoplankton bloom. Furthermore, the model depicts variations in migration depth at fine spatial scales (< 100 km), consistent with those observed across a front by Powell and Ohman (2015), with migration depths varying by up to 200 m between points 10–20 km apart.

In the model, these large-scale, fine-scale and seasonal variations in migration depth are largely caused by the shading of chlorophyll contained in the water column, decreasing the depth of the isolume targeted by migrating zooplankton. Chlorophyll shading leads to differences in migration depth of about 60 m on annual average between the subtropical biome (low-productivity or low-chlorophyll shading) and the subpolar biome (high-productivity or high-chlorophyll shading) and up to 150 m during the spring bloom. In contrast, variations in migration depth caused by differences in surface irradiance between biomes are relatively small (< 40 m), although we note that seasonal variations within biomes are significant (60 to 85 m in the subpolar and convective biomes). We find that large-scale differences in the chlorophyll shading effect across biomes can be reinforced or offset by 100 m locally by fine-scale ocean dynamics associated with eddies and fronts that lead to patchy chlorophyll distribution (Gaube et al., 2014). This patchiness in chlorophyll is consistent with prior studies that showed that fine-scale dynamics could create strong horizontal gradients in chlorophyll content through horizontal stirring (e.g., Martin, 2003) or through the stimulation of phytoplankton growth by local vertical supply of nutrients and changes in light exposure (e.g., Lévy et al., 2018). Also, contrasts in migration depth caused by differences in surface irradiance between biomes are relatively small (< 40 m) but seasonal variations within biomes are significant, ranging from 60 to 85 m in the subpolar and convective biomes.

These results suggest that the ability of migrating zooplankton to sequester carbon in the interior ocean could be limited by chlorophyll shading at both large and fine scales. Productive chlorophyll-rich regions, such as the subpolar gyre at large scales or phytoplankton-rich filaments at fine scales, provide abundant food for zooplankton and are the regions where migrating zooplankton export the most carbon (Lévy et al., 2018; Mangolte et al., 2023; Kaiser et al., 2021; Aumont et al., 2018; Archibald et al., 2019; Nowicki et al., 2022). However, it is also in these regions that chlorophyll shading reduces migration depth the most. The carbon exported by migrating zooplankton in productive regions could therefore be sequestered for a shorter time, as the shallower the migrating zooplankton excrete or respire carbon, the faster it takes for ocean circulation to bring it to the surface, reducing the time it is potentially sequestered in the ocean interior (Boyd et al., 2019).

The model suggests that large-scale ocean transport can further modulate the migration depth, but this effect is confined to deep convective events during which zooplankton can be transported up to 330 m below their isolume. In nature, this process could potentially influence zooplankton migration in the Irminger, Greenland or Labrador seas, where deep convection events occur (Våge et al., 2009; Piron et al., 2017; Sgubin et al., 2017). Our results also suggest that fine-scale ocean dynamics, re-stratifying the mixed layer and causing patchiness in ocean transport (Mahadevan, 2016; Resplandy et al., 2019), can lead to migration depth differences of up to 250 m between two locations 10–20 km apart in the convective biome. However, these results are probably overestimated for two reasons. First, our idealized wind creates an intense downwelling in the convective biome at the northern boundary of the model, leading to the convection of the entire water column (MLD > 3500 m; Figs. 5a, c and A4c), while mixed-layer depths observed in the Labrador, Irminger and Greenland seas are around 200–500 m on average and up to 1500–2000 m at most (Holte et al., 2017; Gonzalez-Pola et al., 2020). Second, in nature, part of the zooplankton community stops diel vertical migration from entering the diapause under the permanent thermocline (600–1400 m) in winter in the North Atlantic (Jónasdóttir et al., 2015), a process not accounted for in our model. At these depths, zooplankton in the diapause could be less affected by deep convection. Therefore, results in the convective biome should be considered with caution.

Existing biogeochemical models that include zooplankton vertical migration primarily use light to control diel migration patterns (Bianchi et al., 2013b; Aumont et al., 2018; Archibald et al., 2019; Nowicki et al., 2022). In contrast, our model introduces a novel prey-driven redistribution mechanism at night, where zooplankton adjust their vertical positions to optimize feeding opportunities by aligning with depths of higher prey concentration, which more accurately captures the critical role of prey availability in shaping zooplankton behavior (Dini and Carpenter, 1992; Beklioglu et al., 2008). However, the model is limited by the lack of predator-led controls, partly because key zooplankton predators are not explicitly represented. This omission could oversimplify zooplankton migration behavior, especially in regions where predator pressure is significant (Bollens and Frost, 1989). In reality, zooplankton DVM is influenced by a balance between avoiding predators and optimizing feeding, with predator pressure often modulating the depth and timing of migrations (Hays, 2003; Bandara et al., 2021). Future improvements could aim to integrate predator–prey interactions more comprehensively to enhance the realism and applicability of the model across various marine environments (Pinti et al., 2019).

This study highlights the role of chlorophyll in setting the zooplankton migration depth in the North Atlantic. Other studies have underlined the importance of other environmental variables such as moonlight (Last et al., 2016) or ice cover (Petrusevich et al., 2016; Flores et al., 2023), modulating the light landscape and isolume depth in the Arctic regions, cloud shadow in the subpolar North Pacific (Omand et al., 2021) or the oxygen concentration limiting zooplankton migration in the upper margin of oxygen minimum zones (Bianchi et al., 2013a). Although estimates are uncertain, climate models predict that most of these parameters will evolve in response to climate change (Kwiatkowski et al., 2020; Notz and Community, 2020; Vignesh et al., 2020). A decrease in ocean productivity or ice cover could deepen the isolume, and ocean deoxygenation could alter the migration depth in oxygen-limiting regions. Moreover, zooplankton size is one of the main features controlling their migration depth (Ohman and Romagnan, 2016). Warmer and less productive waters tend to support smaller zooplankton communities (Barton et al., 2013; Brun et al., 2016), which could dive to shallower depths. The diversity and complexity of these processes make it difficult to predict the evolution of zooplankton migration depths in the future around the globe. However, COBALTv2-DVM captures all these mechanisms, and its implementation in a global ocean model could provide a way of quantifying their importance.

The patchiness of zooplankton migration depth caused by chlorophyll in spring makes it difficult to observe migration patterns and their consequences in the ocean. The patchiness of phytoplankton, which has been discussed for a long time (Mackas et al., 1985; Abraham, 1998; Martin et al., 2002), is easily observable from satellites (Lehahn et al., 2018; McClain et al., 2022). However, chlorophyll heterogeneity spreads to the migration depth of zooplankton, which is difficult to measure. Currently, migration depth is mainly derived from ADCP profiles (Bianchi et al., 2013b; Riquelme-Bugueño et al., 2020; Cisewski et al., 2021) or net measurements (Steinberg et al., 2008; Conroy et al., 2020). These observations are expensive because they are carried out on ships during oceanographic campaigns. Consequently, they are too sparse to resolve and isolate signals of large- and fine-scale spatial variability as well as the seasonal cycle of migration depth in many regions of the ocean (Bianchi et al., 2013a). However, migration signatures have already been observed by ARGO floats (Fig. S3 in Boyd et al., 2019). Developing algorithms to detect these migrations in ARGO floats, as has already been done to detect subduction events (Llort et al., 2018) or deep chlorophyll maxima (Cornec et al., 2021), could be incorporated into integrative observational approaches (Bandara et al., 2021) and could increase the number of observations available. More observations of zooplankton migration are crucial for better understanding their variability, constraining them in ocean models and assessing their impact on biogeochemical cycles.

Appendix A: Materials and methods

Figure A1 Surface forcing seasonal cycle in the double-gyre model. (a) Zonal wind stress following an idealized sinusoidal profile varying by 10 % in amplitude according to a seasonal cycle. (b) Atmospheric surface temperature, (c) sea level pressure, (d) near-surface specific humidity and (e) precipitation derived from the JRA55-do atmospheric reanalysis (Tsujino et al., 2018), with a 3 h time resolution averaged over the period 1958–2020 and zonally over the North Atlantic Ocean (longitude −65 to −20° E, latitude 20 to 60° N). Surface downwelling (f) longwave and (g) shortwave radiation calculated as the zonal mean of the ERA5 atmospheric reanalysis (Hersbach et al., 2020) averaged over the 1979–2022 period and zonally over a band located in the middle of the North Atlantic (longitude −35 to −33° E, latitude 20 to 60° N). (h) Lithogenic material deposition and (i) iron deposition are calculated from the last 30 years of a GFDL-ESM4 preindustrial control simulation spanning the time from 1850 to 2014 in the same way as the physical forcings (Dunne et al., 2020). Each profile is the monthly average of the variables, and the color corresponds to the month.

Figure A2 Comparison between the double-gyre model and the World Ocean Atlas 2018 (WOA18). Zonal mean of (a, b) nitrate concentration, (c, d) phosphate concentration, (e, f) silica concentration, (g, h) oxygen concentration, (i, j) salinity and (k, l) temperature in the double-gyre model (top row) and WOA18 (bottom row, −60 to −20° E and 20 to 60° N).

Figure A3 Map of the annual mean surface eddy kinetic energy (a) in the double-gyre model and (b) derived from sea surface height observations between 1993 and 2021 (Archiving, Validation and Interpretation of Oceanographic Satellite Data – AVISO). The subtropical, jet and subpolar biomes are delimited by solid black lines (mean surface chlorophyll concentration = 0.15 and 0.35 mg m⁻³), while the subpolar and convective biomes are delimited by a dashed line (mean MLD = 200 m).

Figure A4 Seasonal cycle of the (a, b) surface chlorophyll concentration and (c, d) mixed-layer depth in (a, c) the double-gyre model and (b, d) observations as a function of latitude. The seasonal cycle is calculated as the zonal mean over the whole domain in the model and in a North Atlantic box (longitude −65 to −20° E, latitude 20 to 60° N) in the observations. Observations of mixed-layer depths are derived from a compilation of hydrographic sections between 1941 and 2002 (de Boyer Montegut et al., 2004) and surface chlorophyll from satellite observations between 1997 and 2020 (Sathyendranath et al., 2019).

Figure A5 Fraction of zooplankton migrating in the subtropical (blue) and subpolar (orange) biomes in observations (from the COPEPOD database and compiled by Aumont et al., 2018), the double-gyre model and a series of sensitivity experiments. The black line corresponds to the mean and the bars to 2 standard deviations around the mean. The dark-blue and orange bars correspond to the simulation described in the paper and the light-blue and orange bars to the parameter sensitivity experiments.

Figure A6 Observations of the zooplankton gut clearance rate constant versus temperature, (a) binned for each degree of temperature or (b) raw. Digitized data extracted from Irigoien (1998).

Figure A7 Global distribution of vertical migration velocities of zooplankton measured from global acoustic data by Bianchi and Mislan (2016).

Figure A8 Mean daily cycle of the (a) irradiance and (b) visual predation efficiency vertical distribution in the double-gyre model.

Table A1 Major parameters and associated values used in the physical ocean (MOM6) component of the model.

Appendix B: Results

Figure B1 Decomposition of the mechanisms controlling the migration depth of large zooplankton in the (a) subtropical, (b) jet, (c) subpolar and (d) convective biomes. The decomposition includes the effects of large-scale variations in irradiance (yellow bar), chlorophyll shading (green bar) and vertical ocean transport (red bar), together with their seasonal range of variability (black hatch). The effect of transport is nonsignificant (light-red bar) when it is smaller than the vertical resolution of the model grid.

Figure B2 Seasonal decomposition of the mechanisms controlling the medium zooplankton migration depth in (a) the subtropical biome, (b) the jet, (c) the subpolar biome and (d) the convective biome. The lines are the effects of seasonal variations in irradiance (yellow line), large-scale chlorophyll shading (green line) and large-scale vertical transport (red line) relative to a baseline depth set by the mean annual irradiance (isolume in clear water, decreasing from 520 m in the subpolar gyre to 480 m in the convective region). Large-scale effects are superimposed by small-scale spatial variability in chlorophyll shading (green fill) and vertical transport (red fill) induced by eddies and fronts.

Figure B3 Sensitivity experiment on zooplankton migration depth. Change in migration depth and effect of irradiance, chlorophyll and transport for medium (md) zooplankton and large (lg) zooplankton in response to a 50 % decrease or increase in target isolume, half-saturation distance (K_dz) and migration velocity (w). A positive anomaly indicates a deeper migration depth and a negative anomaly a shallower migration depth.