# Dispersion/dilution enhances phytoplankton blooms in low-nutrient waters

Yoav Lehahn, Ilan Koren, Shlomit Sharoni, Francesco d'Ovidio, Assaf Vardi, Emmanuel Boss

► To cite this version:

Yoav Lehahn, Ilan Koren, Shlomit Sharoni, Francesco d'Ovidio, Assaf Vardi, et al. Dispersión/dilución enhances phytoplankton blooms in low-nutrient waters. Nature Communications, 2017, 8, pp.14868. 10.1038/ncomms14868. hal-01502211

HAL Id: hal-01502211

https://hal.sorbonne-universite.fr/hal-01502211

Submitted on 5 Apr 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dispersion/dilution enhances phytoplankton blooms in low-nutrient waters

Yoav Lehahn<fnref n="1" />, Ilan Koren<fnref n="1" />, Shlomit Sharoni<fnref n="1" />, Francesco d'Ovidio², Assaf Vardi³ & Emmanuel Boss⁴

Spatial characteristics of phytoplankton blooms often reflect the horizontal transport properties of the oceanic turbulent flow in which they are embedded. Classically, bloom response to horizontal stirring is regarded in terms of generation of patchiness following large-scale bloom initiation. Here, using satellite observations from the North Pacific Subtropical Gyre and a simple ecosystem model, we show that the opposite scenario of turbulence dispersing and diluting fine-scale (~1-100 km) nutrient-enriched water patches has the critical effect of regulating the dynamics of nutrients-phytoplankton-zooplankton ecosystems and enhancing accumulation of photosynthetic biomass in low-nutrient oceanic environments. A key factor in determining ecological and biogeochemical consequences of turbulent stirring is the horizontal dilution rate, which depends on the effective eddy diffusivity and surface area of the enriched patches. Implementation of the notion of horizontal dilution rate explains quantitatively plankton response to turbulence and improves our ability to represent ecological and biogeochemical processes in oligotrophic oceans.

<footnote n="1">¹Department of Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot 76100, Israel. ²Sorbonne Universités (UPMC, Université Paris 06)-CNRS-IRD-MNHN, LOCEAN-IPSL Laboratory, 4 Place Jussieu, 75005 Paris, France. ³Department of Plant Sciences, Weizmann Institute of Science, Rehovot 76100, Israel. ⁴School of Marine Sciences, University of Maine, Orono, Maine 04469-5741, USA. Correspondence and requests for materials should be addressed to Y.L. (email: yoav.lehahn@weizmann.ac.il) or to I.K. (email: ilan.koren@weizmann.ac.il).</footnote>
A pproximately half of the global primary production is carried out by oceanic phytoplankton¹. To a first order, this immense production of biological matter is controlled by vertical motions that regulate light and nutrient availability, and stimulate the formation of phytoplankton blooms². The biomass accumulated during these blooms is spatially redistributed by horizontal stirring and mixing³⁻⁶. As shown in a number of model studies, in addition to structuring bloom morphology, horizontal stirring may affect different aspects of bloom ecology and biogeochemistry, including regulation of plankton ecosystem dynamics and modulation of primary production levels⁷⁻⁹.

Schematically, the horizontal redistribution process can be divided into two subsequent phases with contrasting effects on tracer distribution. Over short time scales, tracers are stretched and filamented by stirring, gradients are amplified and variance cascades towards smaller scales. Over longer time scales, small-scale turbulence and, ultimately, diffusion mix the tracer with its surroundings and homogenize its distribution. Classically, studies that consider the effect of stirring on phytoplankton blooms refer to the first phase of the process, emphasizing its role in generating mesoscale and sub-mesoscale (~1–100 km, hereafter referred to as fine scale) patchiness following large-scale bloom initiation⁴,⁵,¹⁰. The opposite scenario of fine-scale phytoplankton patches feeding variance back to scales that are larger than the size of mesoscale eddies (typically of the order of 100 km) through dispersion is also expected¹¹. However, the process and its impact on the planktonic ecosystem have largely been overlooked.

A major reason for the limited knowledge about the effect of upscale transfer of plankton variance is the inherent difficulty to track fine-scale seeding events and to follow their evolution as they are advected and dispersed by the oceanic flow. As a result, observational evidence to the effect of horizontal stirring on localized phytoplankton blooms is limited to the case of artificial iron fertilization experiments³, which are limited in their spatial and temporal coverage. Here we address this observational challenge and explore the impact of horizontal stirring on naturally stimulated fine-scale blooms through a detailed Lagrangian characterization of bloom development using multi-satellite products and ARGO float data. This approach allows continuous tracking of bloom dynamics and provides detailed quantification of small-scale spatio-temporal changes in bloom properties. The observations are interpreted using a simple ecosystem model, which simulates the dynamics of planktonic systems embedded within dispersing fine-scale water patches.

Potential ecological and biogeochemical consequences of tracer patch dispersion are of special interest in the oligotrophic waters of subtropical gyres, where a large part of the nutrient supply to the surface waters comes from episodic fine-scale upwelling events¹²⁻¹⁴. We focus our analysis on a region in the North Pacific Subtropical Gyre, where nutrient injection events at these scales have recently been documented¹⁵,¹⁶ and where re-occurring localized summertime blooms are frequently observed in satellite data¹⁷,¹⁸. Satellite observations of summertime blooms in the North Pacific Subtropical Gyre, together with results from the ecosystem model, show that dispersion/dilution of fine-scale nutrient-enriched water patches may regulate the dynamics of nutrients-phytoplankton-zooplankton ecosystems and enhance accumulation of photosynthetic biomass in low-nutrient oceanic environments.

## Results

**Satellite observations.** We first focus on the evolution of one of the largest localized satellite-observed blooms, which took place in the North Pacific Subtropical Gyre during the summer of

2007 (ref. 18). The strong gradients in chlorophyll concentrations allow unambiguous delineation of bloom boundaries according to the chlorophyll contour around its periphery (Fig. 1a–d and Supplementary Fig. 1).

Detailed diagnostics of satellite images (see Supplementary Movie 1) reveals that the bloom was initiated with the emergence of elongated chlorophyll patches covering an area of ~10⁴ km² (Fig. 1a,b). The patches stretched, folded and merged in a continuous manner, resulting in a coherent bloom covering an area of ~10⁵ km² (Fig. 1d). The observed bloom occurred in late summer, when the water column is well stratified and the surface layer is depleted of nutrients¹⁷. This rules out the possibility that the observed patterns resulted from fine-scale stratification of the water column, which was shown to have an important role in triggering spring blooms in the North Atlantic¹⁹. In contrast, these summertime conditions are consistent with submesoscale upwelling of nutrient-rich waters as a trigger for bloom initiation¹²,¹⁴. This is further supported by the fact that the patches analysed here were initiated across distinct sea surface temperature fronts at the peripheries of mesoscale eddies (Supplementary Fig. 2), which are regions that favour the generation of submesoscale vertical velocities²⁰.

The continuous change in bloom morphology suggests that the transition from fine-scale patches was driven by horizontal stirring, followed by mixing³,⁵. To investigate the effect of stirring quantitatively we conducted a numerical simulation in which synthetic particles are advected by the satellite-derived surface currents. The particles were initiated on day 232 (20 Aug 2007, Fig. 1e) according to the morphology of the bloom as delineated by the 0.13 mg m⁻³ chlorophyll contour (Fig. 1b), and were advected by the satellite-derived surface currents until day 290 (21 October 2007; Fig. 1f), when the bloom reached its maximal areal extent (Fig. 1d). After two months of advection, the synthetic particles were dispersed by the stirring process, resulting in dilution of their density. Furthermore, spatial distribution patterns of the bio-optical (that is, chlorophyll) and numerical tracers were remarkably similar (compare Fig. 1d,f), especially considering the limited resolution of the satellite altimetry data used for advecting the synthetic particles. This similarity strongly supports the hypothesis that bloom morphology was structured primarily by horizontal stirring.

Identification of stirring as the factor structuring bloom morphology allowed the construction of unambiguous Lagrangian time series describing changes in the bloom's biophysical properties from satellite data²¹ (Fig. 2). The time-varying patch surface area ($A_p$) is identified as the area bounded by a contour of equal-chlorophyll surrounding the bloom periphery (for example, the 0.13 mg m⁻³ contour in Fig. 1a–d). As $A_p$ depends on the value of the chlorophyll contour, a sensitivity analysis of its temporal evolution was performed for a range of contour levels (0.11–0.16 mg m⁻³; Fig. 2 and Supplementary Fig. 1).

These satellite observations can be used to estimate important physical quantities related to the dispersion process, and in particular the effective eddy diffusivity ($k_e$). Throughout the patch's life time, its surface area increased linearly at an average rate of 2,400 ± 370 km² d⁻¹ (28,000 ± 4,300 m² s⁻¹). The linear increase in $A_p$ is typical for dispersion occurring on a scale that is larger than the size of the local eddies²²,²³. As the plankton patch is dispersed by the turbulent oceanic flow, the rate of its expansion is a measure of $k_e$. Following Garret²², we formulate the linear relationship between $k_e$ and patch expansion rate as follows:

$$k_e = \frac{1}{8\pi} \frac{dA_p}{dt} \qquad (1)$$
<figure><img src="image_3.png" /><figcaption>**Figure 1** | Spatio-temporal evolution of the 2007 bloom. (a–d) Satellite-derived maps of chlorophyll concentrations showing snapshots from the 3 months bloom evolution. Black polygon delineates bloom boundaries as defined by the 0.13 mg m⁻³ chlorophyll contour. (e,f) Density plot showing concentration of synthetic particles advected by the satellite-derived velocity field at the initial and final stages of the numerical experiment (days 0 and 58 of the numerical experiment, corresponding to year days 232 and 290, respectively). The particles were initiated over the area associated with the bloom’s spatial extension on day 232 (that is, b) and were advected without any re-seeding. The positions of the particles in time is provided explicitly by the numerical integrator. Arrows represent the satellite-derived velocity vectors used for advecting the numerical particles. The remarkable similarity between bloom morphology and particle distribution after the 2-month period of particle advection is noteworthy (compare d and f).</figcaption></figure>

Using this formulation, we estimate $k_e$ to be 1,100 ± 150 m² s⁻¹, which is consistent with previous estimates that are based on *in-situ* drifter data²⁴.

Satellite observations of time-varying $A_p$ were used in conjunction with mixed layer depth (MLD) from ARGO floats to quantify the mixed layer chlorophyll biomass associated with the bloom evolution. During its life-time the patch was associated with MLD of 40 ± 10 m, which is typical for summertime conditions in the North Pacific Subtropical Gyre²⁵. Assuming the mixed layer to be homogeneous with respect to chlorophyll²⁶, total chlorophyll biomass was quantified

by integrating satellite-derived chlorophyll concentrations spatially over $A_p$ and vertically over the MLD. Lagrangian time series of the patch evolution shows that the increase in $A_p$ was associated with an order of magnitude increase in total chlorophyll biomass (Fig. 2b). The continuous increase in total chlorophyll biomass indicates that net growth of phytoplankton was positive during the 3-month bloom life time. In contrast to total chlorophyll biomass, averaged chlorophyll concentrations remained almost constant during most of that period (relative change of ~15% during the last 2 months; Fig. 2c).
<figure><img src="image_4.png" /><figcaption>Figure 2 | Satellite-based Lagrangian time series of the 2007 bloom.</figcaption></figure>

(a) Patch surface area. (b) Chlorophyll biomass integrated spatially over the patch surface area and vertically along the MLD as derived from ARGO floats. (c) Chlorophyll concentrations averaged spatially over the patch surface area. Black solid lines and shaded area denote mean and s.d. values, respectively, for different border-delimiting chlorophyll contours (0.11–0.16 mg m⁻³). Blue lines in a and b represent linear and logistic fits, respectively, for the mean time series. Dashed line in c shows background chlorophyll concentrations.

Multi-annual analysis reveals similar patterns (that is, dispersion of fine-scale chlorophyll patches that is associated with relatively constant chlorophyll concentrations and close-to-constant increase in chlorophyll biomass accumulation) in a number of summertime blooms observed in that region (Fig. 3). Lagrangian satellite-based time series for several regional blooms indicate that the amount of particulate organic carbon (POC) accumulated in such blooms can reach the order of 10⁵ tons (Fig. 3e).

**Ecosystem model.** During its dispersion the nutrient-enriched water patch was mixed and diluted with the oligotrophic waters surrounding it. The fact that the linear increase in Aₚ (Fig. 2a) was associated with small temporal changes in chlorophyll concentration (Fig. 2c) indicates that the losses (that is, reduction in concentration) due to dilution were compensated by biological production of phytoplankton²⁷. The competition between

these two effects is especially intriguing as the biological production process is itself influenced by the dilution processes. In principle, dilution can have both a negative (due to reduction in nutrient availability) and a positive (due to reduction in concentration-dependent mortality agents such as viral infection and grazing²¹,²⁸) impact on the phytoplankton's ability to accumulate.

We explore the effect of the dispersion/dilution process on phytoplankton biomass accumulation through an averaged-field numerical model that simulates the dynamics of a nutrient-phytoplankton-zooplankton ecosystem embedded within a water patch whose volume increases with time²⁹ (for a detailed model description, see the Methods section). At time t=0 of the numerical experiment, the patch receives a pulse of nutrients, which represents an episodic event of fine scale nutrient injection. The nutrients are consumed by the phytoplankton, which are grazed by the zooplankton.

We assume that the observed patch is confined within the constant MLD. Accordingly, as the dispersion/dilution processes is associated with a linear increase in patch surface area, it is translated in the model to a linear increase in the volume of the patch, which is driven by a constant inflow of water that do not contain nutrients, phytoplankton and zooplankton. Such volume increase can be translated to a concentration dilution rate. Specifically, we define the horizontal dilution rate (β) as the rate of relative change in patch volume due to the dispersion/dilution processes. Based on equation (1) β can be formulated as:

$$ \beta(t) = \frac{1}{A_p(t)} 8\pi k_e \quad (2) $$

The dilution affects all sink and source processes. For each time step, we solve for the interplay between the coupled production and removal that are regulated by the dilution, and explore the possible solution space and what can be learned from the states to which the systems converges.

We analyse the impact of horizontal dilution on the plankton ecosystem by comparing total biomass of nutrients, phytoplankton and zooplankton accumulated within the patch after the model reaches a steady state, for different initial dilution rates (β₀). A comparison between total biomass levels after the model reaches a steady state for different initial conditions reveals that the potential of a nutrient-enriched water patch to accumulate phytoplankton biomass depends on the balance between characteristic time scales of the physical and biological processes governing its dynamics. Specifically, we distinguish between three scale-dependent regimes that emerge when comparing model results for different values of the non-dimensional parameter γ, defined as the ratio between β₀ and phytoplankton growth rate (μ) (Fig. 4a).

For a given grazing rate (G), when γ ≫ 1 (that is, β₀ ≫ μ), the rapid decrease in nutrient concentration due to dilution prevents the phytoplankton from exploiting the full capacity of the nutrients reservoir and phytoplankton biomass accumulation is reduced compared with the case without dilution. In its extreme, this regime is analogous to the case of patches that are smaller than the critical size required to sustain a plankton patch, as defined by the KISS model³⁰,³¹. When γ << 1 (that is, β₀ << μ) phytoplankton efficiently exploit the nutrient reservoir. However, at a second stage the high phytoplankton concentration stimulates grazing, so that when the system reaches a steady state much of the phytoplankton biomass is converted into zooplankton.

In the third regime, where physical and biological time scales are comparable (that is, β₀ ~ μ), conditions are favourable for accumulation of phytoplankton biomass. In such cases, dilution is
<figure><img src="image_5.png" /><figcaption>**Figure 3** | Dispersed blooms in the North Pacific Subtropical Gyre. (a–d) Maps of satellite-derived surface chlorophyll showing snapshots from the development of fine-scale dispersed blooms observed in four different years. (e) Satellite-based Lagrangian time series of total accumulation of POC for the blooms shown in this figure and for the 2007 bloom discussed in the text.</figcaption></figure>

slow enough to allow efficient utilization of the nutrient reservoir and at the same time fast enough to regulate phytoplankton concentrations such that losses due to grazing are reduced. As phytoplanktons are not starved for nutrients and are relatively protected from grazing pressure, phytoplankton biomass accumulation is enhanced. When these conditions are met, fine-scale phytoplankton patches can thrive for relatively long periods and, ultimately, form blooms that are spread over large areas.

Further examination of the interplay between physical and biological time scales is done by extracting the value of $β_0$ yielding maximum accumulation of phytoplankton biomass ($β_{0\_opt}$, e.g. black star in Fig. 4a) for model runs with different values of $μ$ and G (Fig. 4b). Comparison between the different model outputs shows that for a given grazing rate $β_{0\_opt}$ is linearly correlated with $μ$. The slope of the regression line between $β_{0\_opt}$ and $μ$ is proportional to the value of G.

## Discussion

Using satellite-based Lagrangian time series and an ecosystem model, we have shown how the interplay between governing timescales of plankton dynamics and horizontal dispersion regulates plankton blooms and enhances accumulation of photosynthetic biomass in low-nutrient environments. Owing to the stimulating effect of the dispersion/dilution process, fine-scale nutrient seeding events may be translated into relatively large (~$10^5$ km$^2$) and long-lived (a few months) phytoplankton blooms.

The governing time scale of the dispersion/dilution process is the horizontal dilution rate, $β$, which we define as the relative change in patch surface area due to horizontal turbulent stirring and mixing. Results from a simple ecosystem model embedded within a variable-volume water patch show that ecosystem response to horizontal stirring is dictated by the balance between the physical time scale $β$ and biological time scales $μ$ and
<figure><img src="image_6.png" /><figcaption>**Figure 4** | Plankton response to horizontal dilution in the ecosystem model. (a) Normalized steady-state biomass accumulation of phytoplankton, zooplankton and nutrients (red, black and blue lines, respectively) for different values of the non-dimensional parameter $\gamma$, defined as the ratio between initial dilution rate and phytoplankton growth rate. Phytoplankton growth rate and grazing rate are $1.6 \text{ d}^{-1}$ and $0.12 \text{ d}^{-1}$, respectively. Black star marks maximum phytoplankton biomass yield (associated with $\beta_{O\_opt}$ for the specific combination of growth rate and grazing rate). (b) Optimal initial dilution rate as a function of phytoplankton growth rate. Colours represent different grazing rates. Solid lines show the corresponding linear fits.</figcaption></figure>

G. Optimum conditions for phytoplankton biomass accumulation are found when $\beta_0 \sim \mu$. Under these conditions, phytoplankton can efficiently exploit the nutrient reservoir while being relatively protected from grazing pressure. Thus, by reducing losses due to concentration-dependent factors, the dispersion/dilution process can prolong the life time of fine scale blooms, which otherwise may be rapidly terminated²¹.

The linkage between physical and biological time scales (Fig. 4b) also suggests that the environmental conditions associated with a given dilution rate may be preferential for specific phytoplankton species. Accordingly, fast-growing phytoplankton are expected to thrive in patches associated with high dilution rate, whereas slow-growing phytoplankton are expected to thrive in patches associated with low dilution rate. This explains, quantitatively, the partitioning of phytoplankton functional types into different turbulence regimes, as represented schematically by the mandala of Margalef³². Accordingly, turbulent waters that are characterized by higher values of $k_e$ (and thus, for a given patch size, higher value of $\beta$) would be preferential for faster growing phytoplankton.

Satellite observations reveal that the five dispersed North Pacific blooms discussed here (Fig. 3) are associated with $\beta_0$ values of $0.4 \pm 0.2 \text{ d}^{-1}$, indicating that they are most likely to be found within the favourable regime of $\beta_0 \sim \mu$ and are thus associated with environmental conditions that allow phytoplankton to thrive.

Episodic injection events are an important source of nutrients to surface waters of oligotrophic subtropical gyres. The dispersion and consequent dilution of the resulting fine-scale phytoplankton blooms is therefore likely to be a ubiquitous phenomenon, with substantial impact on the regional ecology and biogeochemistry. As the structure of planktonic ecosystems substantially affect the flux of organic matter to the deep ocean³³, the enhancement of phytoplankton biomass accumulation by horizontal stirring may have an important impact on the efficiency of the biological carbon pump.

In contrast to the case of dilution due to vertical processes³⁴, the ecological and biogeochemical implications of horizontal dilution have barely been explored. Underestimation of the importance of this mechanism may result from an inherent observational limitation: because of dilution the increase in phytoplankton biomass is generally associated with only a marginal increase in surface concentration, which is the

quantity commonly derived from satellite and in-situ measurements. Furthermore, spatial characteristics of the dispersed fine-scale phytoplankton patches are expected to be similar to those resulting from direct cascade of blooms initiated at large scales³⁻⁵. Therefore, identification and quantification of the dilution process can only be achieved by tracking changes in biophysical properties of the dispersed patches in a Lagrangian manner as was done here.

Lagrangian tracking of time varying chlorophyll patches also allows retrieval of important physical quantities related to the dispersion process. Specifically, our results show good agreement between estimated values of $k_e$ (derived from the observed rate of patch expansion) and values derived from in-situ drifter data²⁴. This indicates that dispersion patterns of chlorophyll patches can be used as a direct measure to this important physical quantity and opens the way to systematic retrieval in different parts of the World Ocean. Although there are various ways to estimate tracer dispersion and effective eddy diffusivity in the ocean²³,³⁵⁻³⁷, following chlorophyll patches in a Lagrangian manner using ocean colour satellite data provides a simple and useful observational-based platform for this matter.

Although the fact that biogeochemical properties of phytoplankton blooms may vary in response to changes in spatial distribution of nutrient-supply sources and in rates of horizontal mixing and stirring has previously been acknowledged⁷⁻⁹, the critical effect on the planktonic ecosystem described here has largely been overlooked. The biological effect of horizontal dilution due to tracer patch dispersion is a sub-grid process in large-scale models that do not resolve fine-scale plankton patches. Our approach provides a recipe for parameterizing this effect, by integrating information on the typical sizes of plankton patches and on effective eddy diffusivity, which are the factors that control the horizontal dilution rate. Given the importance of horizontal dilution, such parameterization will improve our ability to represent ecological processes and biogeochemical cycles in nutrient-limited ocean basins.

## Methods

**Satellite data.** Quantification of bloom dynamics is based primarily on a detailed analysis of ocean colour satellite data. To reduce the surface area masked by clouds we integrate observations from three independent sensors: MODIS-Aqua, MODIS-Terra and MERIS. This approach allows continuous tracking of phytoplankton blooms as they are advected and deformed by the currents and provides detailed quantification
of small-scale spatio-temporal changes in bloom properties. The data set comprises 8-day composites of 4 km Level 3 data obtained from the ocean colour data distribution site (https://oceancolor.gsfc.nasa.gov/). From the ocean-colour satellite data, we extract sea surface concentrations of chlorophyll and POC. Both chlorophyll and POC are calculated using band-ratio algorithms, which are based on blue and green regions of the electromagnetic radiation spectrum. The band-ratio algorithms were tested against field measurements of chlorophyll³⁸ and POC³⁹, showing good agreement between satellite-derived and in-situ measured values.

Lagrangian time series of chlorophyll or POC are extracted by averaging the satellite-derived concentrations over the time-varying patch surface area. Total biomass accumulation is derived by integrating concentrations spatially over the patch surface area and vertically over mixed-layer depth extracted from ARGO float temperature profiles (see below). Lagrangian time series of spatially and vertically integrated POC are used as a measure to the amount of organic carbon accumulated in the course of the observed blooms.

Geostrophic surface currents were obtained from the AVISO database (http://www.aviso.oceanobs.com). For the period of our study, the distributed global product combines altimetric data from three satellite missions. The data are available on a 1/3° × 1/3° Mercator grid, with one data file every 7 days. Trajectories used for Lagrangian analysis were computed by integrating the satellite-derived velocities with a Runge-Kutta scheme of the fourth order with a fixed time step of 6 h. The velocities derived from the altimetric data have been linearly interpolated in space and time. Analysis of the velocity field in a Lagrangian manner allows to simulate the formation of spatial structures that are below the resolution of the satellite altimetry data. This is made possible by the fact that small scale patterns in tracer distribution are often associated with the time varying nature of the velocity field, which is captured by Lagrangian analysis tools. Nevertheless, the limited spatial and temporal resolution of the satellite-derived velocity field, together with lack of information wind drift and ageostrophic velocities prevents perfect simulation of the stirring process, resulting in mismatches between spatial distribution patterns of modeled and observed (for example, chlorophyll) tracers.

**Argo floats.** MLD was diagnosed from Argo floats temperature profiles over a 17° × 9° area (197–214°E/27–36°N) overlapping the 2007 bloom. MLD is defined as the depth at which temperature differs from a surface value by 0.2 °C. Argo floats data were collected and have been made freely available by the International Argo Program and the national programmes that contribute to it (http://www.argo.ucsd.edu and http://argo.jcommops.org). The Argo Program is part of the Global Ocean Observing System.

**Ecosystem model for a dispersed nutrient-enriched water patch.** We simulate the effect of patch dispersion and consequent dilution on plankton biomass accumulation, by implementing a simple nutrient—phytoplankton—zooplankton ecosystem model in a variable-volume reactor model²⁶. At time t = 0, the patch receives a pulse of nutrients. Following this single event of nutrient injection, the volume of the patch increases with time in response to constant inflow of seawater that does not contain nutrients, phytoplankton and zooplankton. Assuming that the nutrient-enriched water patch is restricted to a constant depth, the inflow rate (F) of nutrient-free waters is equivalent to the rate of change in patch surface area, Aₚ. Following Garret²¹, we formulate the time-varying Aₚ as:

$$A_p(t)=A_p(0)+8\pi k_c t \quad (3)$$

where k_c is the effective eddy diffusivity. Accordingly, the inflow rate of nutrient-free water is parametrized as:

$$F=8\pi k_c H_{ml} \quad (4)$$

where H_{ml} is the MLD and is taken to be 40 m according to the typical summertime conditions in the North Pacific Subtropical Gyre²⁵.

Phytoplankton concentration is dependent on nutrient concentration and is represented by a Monod growth function. Zooplankton concentration depends on phytoplankton concentration and follows a Holling type II functional response, which assumes that the consumer is limited by its capacity to process food. At time t = 0, the patch contains only trace quantities of phytoplankton and zooplankton. The nutrients injected at time t = 0 are consumed by the phytoplankton, which are grazed by zooplankton. Concentrations of all three ecosystem components are being constantly diluted in response to the inflow of biology-free water into the patch.

The overall dynamics is described by the following equations:

$$\frac{dV}{dt} = F \quad (5)$$

$$\frac{dN}{dt} = -\mu \frac{N}{N+k_N} P - \frac{F}{V} N \quad (6)$$

$$\frac{dP}{dt} = \mu \frac{N}{N+k_N} P - \frac{F}{V} P - GZ \frac{P}{P+k_P} \quad (7)$$

$$\frac{dZ}{dt} = GZ \frac{P}{P+k_P} - \frac{F}{V} Z - mZ \quad (8)$$

where F is inflow rate of nutrient-free seawater into the patch (m³ d⁻¹), V is the patch volume (m³), N, P and Z are concentrations of nutrients, phytoplankton and zooplankton, respectively (mmol N m⁻³), µ is the phytoplankton maximum growth rate (d⁻¹), k_N is the nutrients half-saturation constant (mmol N m⁻³), G is the maximum grazing rate (d⁻¹) and k_P is the half-saturation constant with respect to phytoplankton concentration (mmol N m⁻³). Total biomass accumulation of nutrients, phytoplankton and zooplankton (Fig. 4) are derived by multiplying equations (6)–(8), respectively, by the time-varying patch volume.

The system converges into a steady-state after ~200 days. We simulated the dynamics of phytoplankton, nutrients and zooplankton biomass, and calculated the total biomass accumulated after the system reached a steady state. The steady-state biomass of the nutrients, phytoplankton and zooplankton was normalized to the total carrying capacity of the system, which is defined as follows:

$$CC=P_0 V_0 + N_0 V_0 + Z_0 V_0 \quad (9)$$

where the subscript 0 denotes concentration and volume at time zero.

We assume that all of the biomass is conserved within the mixed layer and ignore losses due to particle sinking, transport by vertical migrators or transformation to dissolved and particulate organic matter. Phytoplankton mortality is neglected under the assumption that zooplankton do not discriminate between live and dead cells. The model does not include losses due to viral and bacterial lysis, and there is no explicit representation of the effect of temperature and light. As a first approximation, we assume uniform fields averaging of temporal and spatial variations in radiation and temperature.

The model was run using MATLAB and its built-in ode45 solver. Parameter values were similar to those used in recent ocean ecosystem model studies⁴⁰,⁴¹ (Supplementary Table 1). Sensitivity tests showed that the presented results are only weakly affected by the choice of parameters (Supplementary Fig. 3) or whether or not grazers are recycled as nutrients.

**Code availability.** The MATLAB codes used for running the ecosystem model are available for use, and can be downloaded from: https://www.dropbox.com/sh/lvdraj2oixcjrr55/AAA8u6Az31jbHmfQgv28dfWca?dl=0.

**Data availability.** All the data used in this research are freely available and may be downloaded through the links detailed in the Methods section.

## References

1. Field, C. B., Behrenfeld, M. J., Randerson, J. T. & Falkowski, P. Primary production of the biosphere: integrating terrestrial and oceanic components. *Science* **281**, 237–240 (1998).

2. Riley, G. A. Factors controlling phytoplankton populations on Georges Bank. *J. Marine Res.* **6**, 54–73 (1946).

3. Abraham, E. R. et al. Importance of stirring in the development of an iron-fertilized phytoplankton bloom. *Nature* **407**, 727–730 (2000).

4. Martin, A. P. Phytoplankton patchiness: the role of lateral stirring and mixing. *Prog. Oceanogr.* **57**, 125–174 (2003).

5. Lehahn, Y., d'Ovidio, F., Lévy, M. & Heifetz, E. Stirring of the northeast Atlantic spring bloom: a Lagrangian analysis based on multi-satellite data. *J. Geophys. Res.* **112**, C08005 (2007).

6. d'Ovidio, F., De Monte, S., Alvain, S., Danonneau, Y. & Lévy, M. Fluid dynamical niches of phytoplankton types. *Proc. Natl Acad. Sci. USA* **107**, 18366–18370 (2010).

7. Martin, A. P., Richards, K., Bracco, A. & Provenzale, A. Patchy productivity in the open ocean. *Global Biogeochem. Cycles* **16**, 9–1–9–9 (2002).

8. Neufeld, Z., Haynes, P., Garcon, V. & Sudre, J. Ocean fertilization experiments may initiate a large scale phytoplankton bloom. *Geophys. Res. Lett.* **29**, 29–1–29–4 (2002).

9. McKiver, W., Neufeld, Z. & Scheuring, I. Plankton bloom controlled by horizontal stirring. *Nonlin. Processes Geophys.* **16**, 623–630 (2009).

10. Abraham, E. R. The generation of plankton patchiness by turbulent stirring. *Nature* **391**, 577–580 (1998).

11. Mahadevan, A. & Campbell, J. W. Biogeochemical patchiness at the sea surface. *Geophys. Res. Lett.* **29**, 32–1–32–4 (2002).

12. Mahadevan, A. The impact of submesoscale physics on primary productivity of plankton. *Annu. Rev. Mar. Sci.* **8**, 161–184 (2016).

13. McGillicuddy, D. J. Mechanisms of physical-biological-biogeochemical interaction at the oceanic mesoscale. *Annu. Rev. Mar. Sci.* **8**, 125–159 (2016).

14. Lévy, M., Ferrari, R., Franks, P. J. S., Martin, A. P. & Rivière, P. Bringing physics to life at the submesoscale. *Geophys. Res. Lett.* **39**, L14602 (2012).

15. Letelier, R. et al. Role of late winter mesoscale events in the biogeochemical variability of the upper water column of the North Pacific Subtropical Gyre. *J. Geophys. Res.* **105**, 28723–28739 (2000).

16. Johnson, K. S., Riser, S. C. & Karl, D. M. Nitrate supply from deep to near-surface waters of the North Pacific subtropical gyre. *Nature* **465**, 1062–1065 (2010).

17. Dore, J. E., Letelier, R. M., Church, M. J., Lukas, R. & Karl, D. M. Summer phytoplankton blooms in the oligotrophic North Pacific Sub-tropical gyre: historical perspective and recent observations. *Prog. Oceanogr.* **76**, 2–38 (2008).
18. Wilson, C. & Qiu, X. Global distribution of summer chlorophyll blooms in the oligotrophic gyres. *Prog. Oceanogr.* **78**, 107–134 (2008).

19. Mahadevan, A., D'Asaro, E., Lee, C. & Perry, M. J. Eddy-driven stratification initiates North Atlantic Spring phytoplankton blooms. *Science* **337**, 54–58 (2012).

20. Calil, P. H. R. & Richards, K. J. Transient upwelling hot spots in the oligotrophic North Pacific. *J. Geophys. Res.* **115**, C02003 (2010).

21. Lehahn, Y. *et al*. Decoupling physical from biological processes to assess the impact of viruses on a mesoscale algal bloom. *Curr. Biol.* **24**, 2041–2046 (2014).

22. Garrett, C. On the initial streakiness of a dispersing tracer in two- and three-dimensional turbulence. *Dyn. Atmos. Oceans* **7**, 265–277 (1983).

23. Sundermeyer, M. A. & Price, J. F. Lateral mixing and the North Atlantic Tracer Release Experiment: observations and simulations of Lagrangian particles and a passive tracer. *J. Geophys. Res.* **103**, 21481–21497 (1998).

24. Zhurbas, V. & Oh, I. S. Drifter-derived maps of lateral diffusivity in the Pacific and Atlantic oceans in relation to surface circulation patterns. *J. Geophys. Res.* **109**, C05015 (2004).

25. de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: an examination of profile data and a profile based climatology. *J. Geophys. Res.* **109**, C12003 (2004).

26. Zawada, D. G. *et al*. A comparison of hydrographically and optically derived mixed layer depths. *J. Geophys. Res.* **110**, C11001 (2005).

27. Okubo, A. in *Spatial Pattern in Plankton Communities* (ed. J. H., Steele) 21–42 (Plenum, 1978).

28. Landry, M. R., Kirshtein, J. & Constantinou, J. A refined dilution technique for measuring the community grazing impact of microzooplankton, with experimental tests in the central equatorial Pacific. *Mar. Ecol. Prog. Ser.* **120**, 53–63 (1995).

29. Dunn, I. J. & Mor, J. R. Variable volume continuous cultivation. *Biotechnol. Bioeng.* **17**, 1805–1822 (1975).

30. Skellam, J. G. Random dispersal in theoretical populations. *Biometrika* **38**, 196–218 (1951).

31. Kierstead, H. & Slobodkin, L. B. The size of water masses containing plankton blooms. *J. Marine Res* **12**, 141–147 (1953).

32. Margalef, R. Life-forms of phytoplankton as survival alternatives in an unstable environment. *Oceanol. Acta* **1**, 493–509 (1978).

33. Siegel, D. A. *et al*. Global assessment of ocean carbon export by combining satellite observations and food-web models. *Global Biogeochem. Cycles* **28**, 181–196 (2014).

34. Behrenfeld, M. J. Abandoning Sverdrup's critical depth hypothesis on phytoplankton blooms. *Ecology* **91**, 977–989 (2010).

35. LaCasce, J. H. Statistics from Lagrangian observations. *Prog. Ocean* **77**, 1–29 (2008).

36. Ledwell, J. R., Watson, A. J. & Law, C. S. Mixing of a tracer in the pycnocline. *J. Geophys. Res.* **103**, 21499–21529 (1998).

37. Lelong, M. P. & Sundermeyer, M. A. Geostrophic adjustment of an isolated diapycnal mixing event and its implications for small-scale lateral dispersion. *J. Phys. Oceanogr.* **35**, 2352–2367 (2005).

38. Werdell, P. J. & Bailey, S. W. An improved bio-optical data set for ocean color algorithm development and satellite data product validation. *Remote Sens. Environ.* **98**, 122–140 (2005).

39. Stramski, D. *et al*. Relationships between the surface concentration of particulate organic carbon and optical properties in the eastern South Pacific and eastern Atlantic Oceans. *Biogeosciences* **5**, 171–201 (2008).

40. Follows, M. J., Dutkiewicz, S., Grant, S. & Chisholm, S. W. Emergent biogeography of microbial communities in a model ocean. *Science* **315**, 1843–1846 (2007).

41. Dutkiewicz, S., Follows, M. J. & Bragg, J. G. Modeling the coupling of ocean ecology and biogeochemistry. *Global Biogeochem. Cycles* **23**, GB4017 (2009).

## Acknowledgements

We acknowledge support by Scott Jordan and Gina Valdez, as well as from the De Botton Center for Marine Science

## Author contributions

Y.L., I.K. and S.S. developed the Lagrangian framework for tracking the dispersed plankton patches, developed the numerical model and analysed the output. Y.L., I.K., S.S., F.D., A.V. and E.B. interpreted the data and contributed to the writing of the manuscript.

## Additional information

**Supplementary Information** accompanies this paper at http://www.nature.com/naturecommunications

**Competing interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**How to cite this article:** Lehahn, Y. *et al*. Dispersion/dilution enhances phytoplankton blooms in low-nutrient waters. *Nat. Commun.* **8**, 14868 doi: 10.1038/ncomms14868 (2017).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

© The Author(s) 2017