## Abstract

Turbulence and coherent circulation structures, such as submesoscale and mesoscale eddies, convective plumes and Langmuir cells, play a critical role in shaping phytoplankton spatial distribution and population dynamics. We use a framework of advection–reaction–diffusion equations to investigate the effects of turbulent transport on the phytoplankton population growth and its spatial structure in a vertical two-dimensional vortex flow field. In particular, we focus on how turbulent flow velocities and sinking influence phytoplankton growth and biomass aggregation. Our results indicate that conditions in mixing and growth of phytoplankton can drive different vertical spatial structures in the mixed layer, with the depth of the mixed layer being a critical factor to allow coexistence of populations with different sinking speed. With increasing mixed layer depth, positive growth for sinking phytoplankton can be maintained with increasing turbulent flow velocities, allowing the apparently counter-intuitive persistence of fast sinking phytoplankton populations in highly turbulent and deep mixed layers. These dynamics demonstrate the role of considering advective transport within a turbulent vortex and can help to explain observed phytoplankton biomass during winter in the North Atlantic, where the overturn of deep convection has been suggested to play a critical role in phytoplankton survival.

## 1. Introduction

Phytoplankton play a central role in the global carbon cycle as it absorbs dissolved inorganic carbon during photosynthesis, some of which is exported to the ocean's interior [1,2]. Ocean dynamics can drive primary production and carbon export at local and global scales [3,4], by controlling biotic and abiotic factors such as light and nutrient limitation, grazing [5] and viruses [6]. Being ultimately driven by light, primary production depends on the net phytoplankton growth between the sea surface and the bottom of the surface mixed layer, the mixed layer depth (MLD). In a seminal paper, Huisman *et al*. [7] investigated the relation between turbulence diffusion, MLD and phytoplankton growth within a homogeneous mixed surface layer and found that growth conditions are not sufficiently described by the classical critical depth [8], defined as the depth where the vertically integrated production balances the net losses within an actively mixed layer. Using a general water column model, they showed that sinking phytoplankton cells manage to persist at intermediate levels of turbulence irrespective of the critical depth [7,9]. Turbulent mixing in the ocean does affect the vertical displacement of phytoplankton cells and, depending on the sinking or buoyancy properties of these cells, different aggregation layers can be found in the water column [10,11].

At high turbulent levels above a critical threshold [7], the increase in light limitation due to a lower average position of the cells in the water column creates unfavourable conditions for growth. On the other hand, low turbulence levels may not be able to sustain sinking cells, thus limiting growth and survival [12]. This relation depends on turbulent mixing and cell sinking rates [13], the latter of which can vary with cell size, species and physiological state [14,15]. Consequentially, the theoretical basis for these critical threshold values of turbulence has been developed in a one-dimensional framework which includes homogeneous vertical eddy diffusion and cell sinking as mechanisms for particle re-suspension and persistence in the water column [7,12]. Eddy diffusion is a common way of including diffusive and advective processes in a single parameter; however, in the ocean both diffusion and advection by themselves can contribute to the mixing of suspended matter, regulating phytoplankton growth and distribution [16]. Indeed, three-dimensional eddy-like structures of various shapes, sizes and intensities are ubiquitous features in the ocean that set the underlying physics driving spatial heterogeneity and the dispersion of marine plankton [17], hence regulating growth and survival. Several physical mechanisms can shape these structures and provide different properties in terms of sinking and dispersion of the suspended matter. For example, shear forces originating from horizontal density gradients can generate mesoscale (10–200 km) and submesoscale (0.1–10 km) eddies [17,18], which typically provide nutrients to the euphotic zone and regulate particle sinking [19–21]. Baroclinic instability in the inner part of the ocean can drive the formation of convective plumes (<1 km) with localized areas of strong downward currents (downwelling) and larger areas of upward currents (upwelling) [22]. On smaller scales, strong and persistent wind conditions can force Langmuir circulation (0.1–200 m) hence generating vortex structures in the vertical with the distinctive long horizontal bands of inert material aligned with wind [23,24].

Although different in their physical nature, all these processes share a similar spatial structure with alternating zones of upwelling and downwelling and zero net transport within a vortex cell [17,22,24]. While these structures are typically three-dimensional, two-dimensional sections can be used to study these systems, capturing relevant advective and diffusive properties of suspended matter [17,18,25–27].

In the present paper, we analysis the effects of turbulent transport in regulating phytoplankton distribution and growth using a novel two-dimensional framework. We use this approach to demonstrate growth and survival of phytoplankton in a turbulent vortex flows at scales commonly found in aquatic ecosystems. We compare our results to previous hypothesis of phytoplankton survival in light-limited environments [7] and apply the model to study the role of deep convection in promoting phytoplankton growth and persistence in winter in the North Atlantic. In these conditions convective overturn has been suggested to maintain phytoplankton cell with sinking rates of several metres per day (electronic supplementary material, appendix C, table C1) within the deep mixed layer and thus to play an important role in sustaining primary production [28,29].

### 1.1. Theoretical framework

We consider the transport and growth of a given phytoplankton population in a turbulent flow field as described by an advection–reaction–diffusion (ARD) equation:
1.1where *P* is the concentration of phytoplankton cells, *w* the phytoplankton sinking velocity and *g*(*z*) the net growth at depth *z*. *K*_{S} is the sub-grid turbulent diffusion, which represents turbulent effect unresolved by the model grid. *U* is the two-dimensional velocity field along the horizontal and vertical dimensions (*x, z*); hence the gradient operator, , is defined on a vertical slice.

Nutrients are assumed to be non-limiting; therefore, net growth is only limited by light availability:
1.2with *g*_{max} being the maximum growth rate, *k _{i}* the light half saturation constant,

*I*(

*z*) the light level at depth

*z*and

*l*a general loss term, which includes respiration, predation and basal mortality.

Light attenuation with depth is calculated according to
1.3where *I*_{0} is the sea-surface irradiance and *k*_{T} is the light extinction coefficient of the water.

The parameter values are presented in table 1.

Within this modelling framework, we split the quasi-random effects of turbulence into a coherent advection field representing convective processes, and spatially varying diffusion representing smaller scale mixing. Specifically, we impose an advective flow described as a time-invariant idealized synthetic flow field of two counter-rotating vortices with zero net transport on the two-dimensional spatial grid (figure 1). The depth of the vortex (*h _{z}*) defines the depth of the surface mixed layer, hence

*h*equals the MLD. The vertical and horizontal velocities are then given by and respectively, with

_{z}*x*being the horizontal position and the stream function (

*ψ*) defined as 1.4

Here *A* is a parameter (units m^{2} s^{−1}) that controls the intensity of the circulation and *u _{z}* is null at the surface and the bottom of the domain. For equation (1.1) we assume no-flux boundary conditions at the surface, periodical boundary conditions on the sides and open boundary conditions at the bottom of the domain.

Sub-grid turbulent diffusion (*K*_{s}) is simulated using Smagorinsky's turbulence closure term [30], which accounts for the sub-grid effects of advection and shear in the vertical and horizontal dimension [31]. A more detailed description of *K*_{s} is given in electronic supplementary material, appendix A.

*K*_{s} accounts for sub-grid turbulence effects on transport, while the advection component is explicitly resolved (equation (1.4)). This formulation is different from the common expression of eddy diffusivity often applied in one-dimensional approaches, where advection and diffusion are combined into one term. To make our results comparable to studies using a single turbulent diffusivity parameter we estimate the overall turbulent diffusivity (*K*_{T}) by converting the mean advection velocity () of the flow field and the sub-grid diffusivity (*K*_{s}) into one diffusivity term (*K*_{T}) [22]:
1.5

Note that *K*_{T} remains spatially structure with lower values at the sea surface, peak values at around 30% of the depth and the lowest values at the bottom (electronic supplementary material, appendix A, figure A2). In convective systems the mean vertical velocity () can be obtained by the ratio between MLD () and a mixing length time scale (*τ*), i.e. [22], with a typical value of *τ* = 12 h [32].

The numerical discretization of equation (1.1) (for more details see electronic supplementary material, appendix B) allows the approximation of the dynamics using a linear time-invariant transition matrix operator (A): 1.6

The linearized formulation of the transport equation is a convenient discrete form of equation (1.1) that allows describing the dynamics of the system by studying the properties of the transition matrix [33], similar to that of matrix population models [34].

We can therefore perform an eigenvalue (and eigenvector) analysis on the transition matrix A to find the stable state solution of equation (1.1); instead of numerically analysing equation (1.1) by iterating it over time until the stable state is reached. Using the eigenvalues and eigenvectors of matrix A it is hence possible to make simple predictions about the persistence and extinction of the plankton community inhabiting a turbulent structure. The dominate eigenvalue of A, *λ*, is related to the long-term net population growth rate (*r*) and, therefore, can be used to represent the stable state solution of the coupled biophysical system i.e. *r* = ln(*λ*). An eigenvalue above one (*λ* > 1) indicates a positive growth rate or bloom condition, while an eigenvalue below one (*λ* < 1) indicates negative growth rate and, in a time-invariant system, will lead to population extinction.

Within the spatial grid *λ* is dependent on the local contribution of one grid cell (*i*) to another grid cell (*j*). The effect of a change in this contribution is proportional to the relative local growth of the destination grid cell *j*, and to the relative abundance of the origin grid cell (*i*) [35]. The relative distribution of the phytoplankton biomass within the two-dimensional domain is proportional to the right eigenvector of the dominant eigenvalue, *v*_{r}, with . In our model, this distribution is the result of phytoplankton growth and transport processes (advection, diffusion and sinking), which can lead to different spatial structures of the biomass aggregation and the relative contribution of local growth to the overall phytoplankton concentration. The latter is proportional to the left eigenvector of the dominant eigenvalue, *v*_{l}, which provides the significance of the different regions to the overall standing stock, a quantity that we define as *reproductive potential*. Indeed, the solution to , indicates how much of the initial conditions are carried over into the dominant eigenvalue. Note that this method yields left and right eigenvectors of the dominate eigenvalue, even if the latter is below one and hence the net population growth rate (*r*) is negative.

This interpretation of the eigenvalues and eigenvectors is a general property common to all matrix differential equation such as equation (1.6), and is analogous to the analysis of the Leslie matrix [36] in modelling the growth of age structured populations in discrete form [34], where the right eigenvector estimates the stable population age structure, while the left eigenvector indicates the reproductive value of each age class.

This framework, which we apply in the following, allows to explicitly resolve the combined effects of light-dependent growth, cell sinking, turbulent diffusion and advective transport to determine the spatial distribution of biomass and net population growth rates of phytoplankton in the model.

## 2. Results

### 2.1. Spatial structure in a vortex field

Assuming zero advection velocity ( = 0), the dynamics are similar to a one-dimensional model, hence the spatial distribution of phytoplankton is horizontally homogeneous (figure 2*a,b*). The highest aggregations occur at the surface when the cells are neutrally buoyant (*w* = 0 m d^{−1}, figure 2*a*) and at the bottom at higher sinking rates (*w* = 10 m d^{−1}, figure 2*b*).

At > 0, neutrally buoyant cells accumulate both in upwelling and downwelling areas (figure 2*c*), with the latter becoming more elongated at increasing flow velocity (figure 2*e,g*). Sinking cells aggregate in the centre of the vortices (figure 2*d*), while further increases in flow-field velocity lead to a more homogeneous distribution throughout the domain with elevated concentrations in the upwelling areas (figure 2*f,h*). Numerical simulations of equation (1.6) confirm those results (see electronic supplementary material, appendix B, figure B1).

In all cases with a non-zero velocity field, the *reproductive potential* (left eigenvector, *v*_{l}) is maximized in the regions of local upwelling independently of the sinking velocity (figure 3*c–h*). However, sinking phytoplankton exhibit higher values of *v*_{l} in the centre of the vortex (figure 3*d,f,h*), whereas neutrally buoyancy phytoplankton shows the lowest values in the centre (figure 3*c,e,g*).

These spatial dynamics rely on the existence of the turbulent vortex. When parametrizing the advective flow into a mean vertical profile of *K*_{T} (electronic supplementary material, appendix A, figure A2) and applying it homogeneously in the horizontal dimension, the vertical distribution of the phytoplankton concentration changes (electronic supplementary material, appendix A, figure A3) and as expected, phytoplankton concentrations are homogeneous in the horizontal dimension.

### 2.2. Vortex effects on sinking plankton cells

The strength of the velocity field has important implications for the maintenance of phytoplankton cells within the actively mixed layer. Both MLD and advection velocity in the vortex contribute to shaping the range of phytoplankton sinking speeds allowing for positive population growth (figure 4).

Generally slower sinking cells have higher population growth than faster sinking cells. When advective transport is weak, lower sinking velocity leads to higher growth rates relating almost linearly to the MLD (figure 4*a*).

With increasing advection, the negative effect of increasing light extinction on the net population growth rate decreases (electronic supplementary material, appendix A, A5). Interestingly, the magnitude of decrease in growth rate shifts with increasing advective strength from shallower to deeper MLDs (figure 4*b–h*). With increasing fluid velocity the role of sinking in regulating growth also becomes weaker allowing cells with higher sinking rate to maintain positive growth rates (figure 4*b–h*). However, at high levels of advection the growth rate again scales linear with sinking velocity and the MLD (figure 4*i*).

These dynamics can be explained by the cells tendency to aggregate in regions of low velocity, which are those in the inner part of the vortex cells (figure 2*d*).

### 2.3. Critical conditions for phytoplankton growth

By calculating a single parameter of turbulent diffusivity (*K*_{T}) from the advective flow velocities and sub-grid turbulence processes (*K*_{s}), our results become comparable to the results by Huisman *et al*. [7,9]. Under these conditions the model reproduces boundaries for phytoplankton growth that are consistent to the results for a one-dimensional water column [9] (figure 5*a*).

The highest asymptotic growth rates for the phytoplankton population (*r* = 0.43 d^{−1}) is found in shallow mixing areas (∼*z* = 10 m) at levels of intermediate and high diffusivity (*K*_{T} = 10^{−2}–10^{0} m^{2} s^{−1}). At intermediate turbulent diffusion and for deeper mixing layers, *r* remains constant up to the maximum depth tested (1000 m) (figure 5*a*). Hence, when using *K*_{T}, phytoplankton growth above a certain MLD only scales with diffusivity and becomes independent of the MLD.

However, when we apply an explicit turbulent flow-field and sub-grid diffusivity *K*_{s} in two dimensions, some of the conditions for phytoplankton growth are changed (figure 5*b*). For growth rates above a certain threshold of turbulence, the scaling is now dependent on both, the MLD and mixing generated by advection and diffusivity. Thus, contrary to previous findings where positive growth rates were limited to intermediate turbulence levels, positive growth rates can now be maintained at high mixing depth when turbulence is also high (figure 5*b*).

### 2.4. Growth in deep convection systems

To study phytoplankton growth in deep convection systems we run the model using a constant mixing-length time scale, *τ* (see model description), with the related turbulent energy dissipation rates comparing well with observations (electronic supplementary material, appendix A, figure A1). We analyse phytoplankton growth rate at different cell sinking rates and MLDs. The overall growth rates in deep, highly turbulent waters are low with the highest growth rates being simulated for slow sinking cells in shallow mixed layers (figure 6). However, because mean available light decreases with increased convective mixing depth and increasing advective transport with increasing depth compensates for cell sinking, the balance between these two dynamics leads to an optimal depth where the highest sinking rates can be sustained. Using the parameter values in table 1, this depth was found around 300 m, with a maximum sinking speed at about *w* ≈ 6 m d^{−1} sustaining a positive growth of the phytoplankton population (figure 6). When compared with observations, the range of positive growth predicted by the model fits well with measured sinking rates of species found in convective regimes in the North Atlantic (electronic supplementary material, appendix C, table C1).

## 3. Discussion

Turbulent structures in the ocean play an important role in driving population dynamics of phytoplankton and their spatial structures [3]. They do so not only through regulating nutrient availability [19] but also by directly impacting on their spatial distribution [16] and hence their access to light. Our results confirm that the interplay between growth, sinking and turbulent transport can drive phytoplankton blooms and affect phytoplankton spatial distribution within the mixed layer. Moreover, explicit consideration of the turbulent vortex flow has a significant effect on the net phytoplankton growth and survival. We apply a numerical technique based on the analyses of a transition matrix that is commonly used in demographic studies in ecology [34]. We extended the method to analyse spatial dynamics on a two-dimensional grid and applied it to study the spatially structuring effects of vortices of turbulent flow. In addition to the net population growth rate, the method provides quantitative information on biomass accumulation (figure 2) and growth potential (figure 3) of phytoplankton in the spatial domain, using the right and left eigenvectors, respectively.

Using this method, this study demonstrates the spatially structuring effect of turbulent vortex flow on phytoplankton growth, and biomass accumulation can explain the existence of a viable phytoplankton stock with significant sinking rates commonly found in the North Atlantic convective systems during winter [28].

### 3.1. Spatial structure

Turbulent vortices in the ocean can be relatively coherent structures. Convective plumes, Langmuir circulation, and meso- and submesoscale vortices can continue to exist for hours [37], days [38] or months [17,39], respectively. These timescales are considerably longer than the simulated values of the Kolmogorov time scales (*τ _{η}*), defined as , which depend on the velocities and depth of the vortex and are on the order of seconds or minutes in these systems.

In convective systems the observed turbulent energy dissipation rate (*ɛ*) ranges on the order of 10^{−10} to 10^{−6} W kg^{−1} [37,40–42]. A similar range was found by Plueddemann *et al.* [43], who measured *ɛ* in Langmuir circulation finding values ranging in between 10^{−10} and 10^{−7} W kg^{−1}, while large eddy simulations of Langmuir circulation [44] suggested *ɛ* to range between 10^{−8} and 10^{−6} W kg^{−1}. Knowledge about *ɛ* in large mesoscale eddies is limited [45]. However, for submesoscale eddies observations in the Southern Ocean Garabato and co-authors [46] suggested values around 10^{−10} W kg^{−1} and model simulation [47] suggested values of *ɛ* reaching up to 10^{−8} W kg^{−1}. Using the mixing length time scale *τ* = 12 h in our model, we obtain similar values of *ɛ* ranging from approximately 1.25 × 10^{−10} W kg^{−1} at *h _{z}* = 100 m to approximately 5 × 10

^{−08}W kg

^{−1}

*h*= 2000 m (electronic supplementary material, appendix A, figure A1).

_{z}For sinking particles several spatial distributions in relation to vertical vortex cell have been suggested.

In our simulations the area of accumulation of sinking cells changes with the strength of the flow field.

At lower flow-field velocities cells accumulate in the centre of the vortex (figure 2*c*), supporting the idea of a retention zone, where sinking particles are being trapped, hence inhibiting sedimentation [25].

More recent work suggests that such trapping is only possible if more complex flow structures are considered [48]. Nevertheless, even if sedimentation of cells occurs continuously, the persistence of a viable population ultimately depends on relation of phytoplankton growth and sedimentation [7], an argument that is supported by model studies [18,49].

With increasing flow-field velocities cells are first been spread out more homogeneously (figure 2*f*) and eventually accumulate in the upwelling areas of the vortex (figure 2*h*). Interestingly, both of these accumulation patterns also have been suggested previously. Theoretical considerations suggest that sinking cells mainly accumulate in upwelling regions [50], whereas analysis of individual trajectories of slowly sinking particles in a Lagrangian model showed a random distribution of particles within the mixed layer [27].

Based on observations [23,50] and numerical modelling [27] vortex-like turbulent flows have been suggested to induce accumulation in the downwelling regions of Langmuir cells in neutrally and positively buoyant suspended cells.

Our results are consistent with these findings and suggest that neutrally buoyant cells do indeed accumulate in the downwelling areas of simple vortex structures (figure 2); however, we also found accumulation of cells in the upwelling area. This pattern changed very little with increasing flow velocity.

Both for neutrally buoyant and sinking cells our model suggests that generally the maximum values of the reproductive potential, the relative contribution of the local growth to the overall phytoplankton concentration, are located close to the surface in the upwelling areas (figure 3). In these areas cells will subsequently be transported horizontally near the surface, before again being transported to greater depth, hence maximizing the amount of light obtained.

These results show that the different patterns of accumulation of sinking cells within a turbulent vortex structure can be explained by different turbulent flow velocities. We observe a transition from accumulation in the centre of the cell towards an accumulation in the upwelling areas, with a transitional more homogeneous distribution. The dependency of these spatial distributions is confirmed when comparing these findings with simulations using the same mean vertical *K*_{T} (electronic supplementary material, appendix A, figure A2), but a homogeneous diffusivity profile in the horizontal dimension (electronic supplementary material, appendix A, figure A3). Under these conditions phytoplankton concentrations are homogeneous in the horizontal dimension and the highest concentrations are either found at the surface or at the bottom of the domain. It is worth noting, however, that we only considered passively sinking cells at a constant rate, while mobile phytoplankton, in particular gyrotaxis [11], as well as variable sinking rates [49], have been shown to influence spatial distribution and intensive patchiness. In fact, phytoplankton cells have been shown to react to turbulent cues by altering their swimming direction in order to avoid exposure to highly turbulent layers [51]. Nevertheless, our model indicates that increased phytoplankton concentrations can occur along the horizontal dimension (figure 2), suggesting that the described dynamics could potentially play a role in the observed fine vertical phytoplankton distribution [10,52].

### 3.2. Bloom conditions

In our model the net population growth rate generally decreases with increasing advection (figure 4). Experimental results on the influence of mixing on the net growth rate are ambiguous [53,54], suggesting that multiple factors play an important role in shifting between a positive and negative feedback. Ruiz *et al.* [55] found that increased mixing leads to increased settling of cells, hence reducing net growth. Our results indicate that for a given sinking velocity the net population growth depends on the combination of turbulence (advection and diffusion) and MLD (figure 5*b*).

If turbulence exceeds a critical threshold, turbulent diffusion will dominate sinking and growth yielding a negative net population growth. Similarly, if turbulence falls below a minimal value, sinking cells are not sustained in the water column thus not supporting a viable phytoplankton population.

This general picture of a window of turbulence for phytoplankton blooming conditions, as described by Huisman *et al*. [7,9], is confirmed when only considering homogeneous turbulent diffusivity (*K*_{T}) (figure 5*a*), and our analysis reproduces previous results, where beyond a certain depth the region of positive growth only depends on the level of mixing [7,9,12].

Under these conditions, the maximum critical depth can be found at around 55 m for values of *K*_{T} exceeding ∼40 cm^{2} s^{−1} (figure 5*a*). While this value is significantly higher than the background diffusivity commonly found in the inertia of the ocean [56], submesoscale eddies [57], convective plumes [22,58] and Langmuir circulation [59] are highly turbulent regimes, where *K*_{T} frequently reaches levels that are orders of magnitude higher. These systems are characterized by a similar vortex cell, with comparatively high advective transport within the cell.

Indeed, when an explicit representation of a two-dimensional vortex structure is introduced we identify important changes in phytoplankton spatial distribution and survival. Above a certain level of turbulence, sinking phytoplankton accumulates in the centre of the vortices, allowing it a higher position within the water column, and hence reducing light limitation in deep waters (figure 2). This allows deeper waters to sustain not only higher rates of cell sinking, but also higher net growth rates (figure 4). These conditions suggest that the region of positive growth in deeper water does not only depend on the level of mixing, but also on the MLD. Hence, growth is driven by a combination of mixing intensity and mixing depth. In fact, net population growth scales linear with the ratio of mixing and water column depth (figure 5*b*), which is explained by the impact that the flow field has on the spatial structure of phytoplankton distribution.

### 3.3. Relevance for deep convection systems

It has been suggested that highly turbulent and deep mixed layers are not able to sustain positive phytoplankton growth [9]. However, observations of viable phytoplankton stocks in such waters [28,29,60] challenge those theoretical findings. Backhaus *et al*. [28] suggested that this discrepancy stems from the missing representation of convective overturn, when using only turbulent diffusion to parametrize turbulence. Similarly, upward convective transport has been suggested to sustain sinking diatom cells during winter in the Irminger Sea [29].

From theoretical consideration of deep convective regimes, one would expect *K*_{T} to increase with increasing depth as overturning is driven by a buoyancy loss at the surface, resulting in a constant mixing length time scale (*τ*). Indeed, scaling of *K*_{T} with mixing depth (*h _{z}*) is supported by observations of turbulent diffusivity from the Mediterranean, Greenland and Labrador Seas and by theoretical scaling of turbulence with negative net surface heat flux (electronic supplementary material, appendix A, figure A4).

Both, for the observations and the theoretical scaling, *K*_{T} increases with increasing depth, exceeding the upper turbulence limited at which a positive phytoplankton growth rate has been suggested previously (figure 5*a*). However, when considering advection explicitly the scaling of phytoplankton growth with turbulence (figure 5*b*) suggests a positive growth rate at the observed levels of turbulence, thus helping to explain the observed phytoplankton biomass during winter in such systems.

Our results further show that when advection is considered explicitly, phytoplankton cells with significant sinking rates can be sustained within a deep highly turbulent mixed layer and form a growing phytoplankton community (figure 6).

Hence, our simulations lend support to the ‘phyto-convection’ hypothesis [28], which states that deep convection during winter sustains a viable phytoplankton community by superimposing cell sinking and frequently returning cells into the euphotic zone and hence providing the inoculum for the following spring bloom [61].

### 3.4. Diversity

Spatial heterogeneity, or patchiness, generated by turbulent vortices in the ocean can increase diversity [62]. Our model is not able to address the impact of spatial heterogeneity on the diversity of phytoplankton directly, as the model is constructed for non-limiting nutrient conditions (typical of deep convection regions and pre-bloom phases) and does not include self-shading. Therefore, different phytoplankton populations will not be competing in the model and results of ‘multi-trait’ simulations with a set of different phytoplankton groups (e.g. in their sinking rates) will produce results equivalent at steady state.

Nevertheless, the model provides indications that in turbulent vortex structures phytoplankton cells of different sinking speeds can be spatial separated. This suggests that the interplay between cell sinking and the turbulent structures can lead to different heterogeneous phytoplankton distributions in the water column. Moreover, we also show that above a certain level of turbulence, the range of sinking velocities that can produce a viable phytoplankton population is relatively large (figure 4); hence suggesting that species with different size and shapes could all produce a growing population.

### 3.5. Advantages and limitations of the modelling framework

The modelling framework presented in this study provides steady-state solutions for an ARD equation (equation (1.1)). In simple spatial explicit biophysical models, e.g. NPZD-type water column models, such steady-state solutions are commonly found by numerically iterating through time until a steady stage is reached. In comparison, the method presented here provides the same result (see electronic supplementary material, appendix B, figure B1) using an eigenvalue analysis, hence an analytical approach which does not only ensures numerical stability but also reduces computing time. Huisman *et al*. [12] in their appendix presented a fast algorithm providing quantitate binary (bloom or non-bloom) information regarding blooming conditions. Compared to their numerical recipe, our approach yields not only actual values of the net population growth rate, but additionally provides the spatial distribution and the *reproductive potential* of phytoplankton, while maintaining faster computing time.

Providing only stable state solutions, our approach neglects temporal dependencies, making it unable to resolve transient dynamics, time-dependent feedbacks (e.g. nutrient limitation), and temporally varying parameters such as variable cells sinking. Similar to most water column models using relatively simple ARD equations, our approach does not account for feedbacks of cell sinking on turbulent flow, which can be accounted for by modelling the Navier–Stokes equations explicitly.

Prescribing fixed values for the biological parameters, the model does not capture the cells' adaptive capabilities, which can play an important role in the survival and seasonal succession of phytoplankton [63]. For example, by applying a fixed value for *w* in equation (1.1), the model ignores the ability of actively controlling the swimming direction in reaction to changing turbulent conditions [51] and to move in between the surface layer and the deeper ocean, as done during dial vertical migration [64]. Such behaviour allows cells to obtain nutrients from below the nutricline, optimize energy allocation and reduce top-down control. Not being able to capture these dynamics our modelling results can, therefore, be considered to produce conservative estimates, as cell sinking out of the surface layer are simply lost to the standing stock.

Further, this stationary approach does not account for temporal variability in mixing strength, which can affect phytoplankton growth [3,17] and may induce blooms in the absence of stratification [65].

However, given that a reduction in mixing leads to improved growth conditions, the assumption of a stationary vortex is likely to underestimate the average growth condition and, therefore, provides a conservative estimate. Moreover, our model is strictly valid in light-limited environments where nutrients are not limiting. To solve such nutrient-dependent dynamics a system of differential equations would be needed and the numerical procedure we use to obtain the transition matrix will not be applicable any longer.

Nevertheless, our results can for example shed light on the winter productivity of phytoplankton during the convective mixing period in the North Atlantic (figures 5 and 6). In such system, prior to the onset of thermal stratification, growth is mainly light-limited and phytoplankton concentrations are low so that self-shading is likely to be of lesser importance. Under these conditions the potential growth can be well approximated by light availability and the dynamics approximated by the simple advection diffusion reaction used in our study.

## Data accessibility

Additional data and analyses are made available in the electronic supplementary material, appendix A, B and C.

## Authors' contributions

C.L. designed the study, drafted the manuscript, wrote the model code and performed model analysis; A.W.V. helped with the modelling and commented on the manuscript; P.M. designed the study, contributed to the writing of the manuscript and helped with the modelling.

## Competing interests

We declare we have no competing interests.

## Funding

C.L. was partially financially supported by the FP7 program EURO-BASIN (grant no. 264933). The VKR Centre for Ocean Life (P.M., A.V.) is supported by the Villum foundation.

## Acknowledgements

We thank the referees for their valuable comments and suggestions, Dag L. Aksnes for useful comments on the manuscript and Christoffer Klausmeier, Jorn Bruggeman and Adrian Martin for helpful comments on previous versions of the manuscript.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3911878.

- Received June 19, 2017.
- Accepted October 9, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.