The effect of intermittent upwelling events on plankton blooms

TL;DR

This study investigates how intermittent nutrient upwelling events influence plankton blooms, revealing that the hydrodynamic flow structure critically determines whether nutrients promote harmful algal blooms or are washed away.

Contribution

It demonstrates that the response of plankton to nutrient pulses depends on mesoscale flow dynamics, explaining the variability in bloom formation and observational inconsistencies.

Findings

Nutrient upwelling can either trigger blooms or be washed out.

Flow vortices determine nutrient retention or dispersal.

Predicting HABs requires understanding flow-nutrient interactions.

Abstract

In the marine environment biological processes are strongly affected by oceanic currents, particularly by eddies (vortices) formed by the hydrodynamic flow field. Employing a kinematic flow field coupled to a population dynamical model for plankton growth, we study the impact of an intermittent upwelling of nutrients on triggering harmful algal blooms (HABs). Though it is widely believed that additional nutrients boost the formation of HABs or algal blooms in general, we show that the response of the plankton to nutrient plumes depends crucially on the mesoscale hydrodynamic flow structure. In general nutrients can either be quickly washed out from the observation area, or can be captured by the vortices in the flow. The occurrence of either scenario depends on the relation between the time scales of the vortex formation and nutrient upwelling as well as the time instants at which upwelling pulse occurs and how long do they last. We show that these two scenarios result in very different responses in plankton dynamics which makes it very difficult to predict, whether nutrient upwelling will lead to a HAB or not. This explains, why observational data are sometimes inconclusive establishing a correlation between upwelling events and plankton blooms.

Tables

Table 1. Table 1: The steady state values in gC m -3 for the systems subjected to k d subscript 𝑘 𝑑 k_{d} .

	Sys.(I)	Sys.(II)
$k_d$	$0.015$	$0.0045$
$N$	$0.02$	$0.034$
$P_N$	$0.04$	$0.007$
$P_T$	$0.024$	$0.006$
$Z$	$0.028$	$0.016$

Table 2. Table 2: Spatial averages in gC m -3 for the system without upwelling.

	Sys.(I)	Sys.(II)
	$0.025$	$0.0286$
	$0.042$	$0.0084$
	$0.019$	$0.0010$
	$0.005$	$0.0054$

Table 3. Table 3: Parameters used in the hydrodynamic flow model (for details of parameters see   [ 29 ] ).

Parameter	Symbol	Used value
Island radius	$r$	$25$ km
Horizontal main flow velocity	$u_0$	$0.18$ m s-1
Velocity of the Ekman flow	$u_E$	$0.018$ m s-1
Vortex strength	$\omega$	$55\cdot {10}^3$ km2 s-1

Table 4. Table 4: Functional responses used in Eq.(1) of the main text

Nutrient uptake by non-toxic species of phytoplankton	$f_N\left(N\right)$	$\frac{N}{e_N+N}$
Nutrient uptake by toxic species of phytoplankton	$f_T\left(N\right)$	$\frac{N}{e_T+N}$
Growth rate limitation due to light attenuation	$g(P_N,P_T)$	$\frac{a}{b+c{P}_N+c{P}_T}$
Feeding rates of the Zooplankton	$h(P_N,P_T)$	$\frac{1}{\left({\mu }^2+P_N^2+P_T^2\right)}$

Table 5. Table 5: The values used are taken from ranges given in  [ 47 ]

Parameter	Symb.	Sys.(I)	Sys.(II)	Units
a/b maximum daily nutrient uptake	$a$	$0.2$	$0.2$	m-1day-1
Light attenuation by water	$b$	$0.1$	$0.1$	m-1
Phytoplankton self-shading coefficient	$c$	$0.4$	$0.4$	m2 gC-1
Mortality rate of Zooplankton	$d$	$0.065$	$0.065$	day-1
Half-saturation rate for $N$ uptake of $P_N$	$e_N$	$0.02$	$0.02$	gCm-3
Half-saturation rate for $N$ uptake of $P_T$	$e_T$	$0.1$	$0.1$	gCm-3
Respiration rate of $P_N$	$r_N$	$0.1$	$0.1$	day-3
Respiration rate of $P_T$	$r_T$	$0.05$	$0.1$	day-3
Conversion rate of nutrients into $P_N$	${\theta }_N$	$0.4$	$0.4$
Conversion rate of nutrients into $P_T$	${\theta }_T$	$0.8$	$0.4$
Phytoplankton sinking rate	$s$	$0.08$	$0.08$	day-1
Growth efficiency of $Z$ due to $P_N$	${\alpha }_N$	$0.25$	$0.5$
Growth efficiency of $Z$ due to $P_T$	${\alpha }_T$	$0.2$	$0.2$
$Z$ excretion fraction	$\beta$	$0.33$	$0.33$
Excretion factor of $Z$	$\gamma$	$0.5$	$0.5$
Maximum grazing rate of $Z$	$\lambda$	$0.65$	$1.3$	day-1
Grazing of $Z$ half saturation constant	$\mu$	$0.02$	$0.02$	gC m-3
Intensity of grazing on $P_T$	$\varphi$	$0.5$	$0.05$

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}t}$$

$$+\beta h(P_N,P_T)\left[\lambda (1-\varphi )P_N^2+\lambda \varphi P_T^2\right]Z+\gamma dZ,$$

$$\frac{\mathrm{d}{P}_N}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{P}_T}{\mathrm{d}t}$$

$$\frac{\mathrm{d}Z}{\mathrm{d}t}$$

$$\frac{\partial C}{\partial t}+\vec{u}\nabla C=F_c+D{\nabla }^{2}C,\text{with}C\in [N,P_N,P_T,Z],$$

$$\begin{array}{rl}\dot{x_1}& =x_2,\\ \dot{x_2}& =-x_1^3-2x_1{x}_3+x_1{x}_5-\mu x_2,\\ \dot{x_3}& =x_4,\\ \dot{x_4}& =-x_3^3-{\nu }_1{x}_1^2+x_3{x}_3-{\nu }_2{x}_4,\\ \dot{x_5}& =-{\nu }_3{x}_5-{\nu }_4{x}_1^2-{\nu }_5(x_3^2-1),\end{array}$$

$$k(t^{\prime })=(k_{\text{up}}-k_d)|x_1(t^{\prime })|+k_d.$$

Full text

Keywords: upwelling, eddies, harmful algal blooms.

1 Introduction

Coastal regions susceptible to harmful algal bloom (HAB) events are often subjected to upwelling [1]. Due to this upwelling nutrient-rich deep waters are transported into the euphotic zone and this inflow fosters favorable conditions for the growth of algae [2]. As recent studies notice a significant increase of the number of harmful algal bloom events in the whole world [3, 4, 5], it becomes imperative to understand the interplay between the biotic and physical factors that work as their trigger.

Lateral mixing and stirring by the hydrodynamic flow redistributes the nutrients and the suspended microorganisms, shaping the spatial heterogeneity of the marine ecosystem at different scales [6], leading to plankton blooms, which exhibit a non-uniform distribution in space, referred to as “patchines”. This non-uniformity was ubiquitously detected around the globe by satellite imagery [7, 8, 9, 10] and by samples along ship transects [11, 12, 13]. In a seminal work by Abraham [14], a very simple model of turbulent transport was able to reproduce this spatial heterogeneity in the plankton distribution and its statistical properties, such as spectra. This model shows that advection by ocean flows on the mesoscale (10 - 100 km) can spatially distort the concentration of plankton leading to the development of small spatial patterns and thin filaments. Subsequent theoretical studies have observed that the ratio between biological and hydrodynamic flow time scales has a non-trivial impact on how plankton is distributed spatially [15, 16]. The flow field influences not only the spatial distribution but also the abundance of plankton. Further studies have shown that it is possible to trigger or suppress HABs by tuning the flow to the biological timescales [17, 18]. Therefore the hydrodynamics plays a central role for phytoplankton ecosystems, not only with respect to its spatial patterns but also to the inter-specific interactions, establishing so called “fluid dynamical niches”, which provide particular growth conditions for certain species [19, 20]. Coherent structures of the flow field, such as, for instance eddies play an important role influencing the biological processes in the ocean. The recent advances in detecting and tracking eddies in the ocean have shown that they often are long lived. Notably they can trap fluid and the whole community of plankton and bacteria inside, which affects the diversity and dominance structure of phytoplankton species observed in the system [21, 20]. On the one hand, the species in the almost isolated ecosystem inside the eddy are subjected to competitive pressure. On the other hand, theoretical models speculate that due to this trapping the organisms can also be sheltered inside the eddy from predators or competitors, a mechanism proposed as a possible explanation for the coexistence of species [22, 23, 24, 25]. Both effects are a direct result of transport barriers established by the flow field.

The productivity enhancing effect of coastal upwelling is also shown to be strongly affected by the presence of coherent structures in the flow field [26, 27, 28]. The eddies mix and disperse nutrients, while also taking them away from the coastal region. This leads to a decrease in the primary production near the shore, as was recorded for eastern boundary upwelling systems [28]. Furthermore, these nutrients while being transported offshore by the eddies may also trigger the growth of the associated phytoplankton. Theoretical works have shown that eddies in this case work as incubators for growth by sustaining favorable environmental conditions. These models emphasize the importance of the role of biological and hydrodynamic timescales in triggering plankton blooms and specifically HABs in this scenario [29, 30, 31, 32]. However, these studies, have so far only analysed the conditions of an upwelling which is constant in time. Nevertheless, upwelling itself is not a steady process since it depends on winds and seasonality, being therefore highly intermittent. Furthermore, from observations of HABs in nature, it is not always possible to correlate the strength and the duration of an upwelling event and the occurrence and magnitude of HABs. While it was possible to establish such a direct relation for some species (e.g. diatoms of genus Pseudo-nitzschia [33, 34] and some dinoflagelates  [35]), for others a more complex chain of events appears to be driving the outcome [36]. Moreover, the major challenge consists in finding out how the occurrences of HABs and the episodic upwelling events are associated on a local scale [36, 37].

In this work we analyse the impact of intermittent upwelling events on phytoplankton growth and changes in dominance patterns in the presence of mesoscale hydrodynamic structures for a biological system with three trophic levels. We modify the reaction-advection-diffusion model introduced by Sandulescu et al. (2006) for the area around the Canary islands  [29], that couples advection by a vortex street behind an island with a model of plankton dynamics. As in  [29] we also choose to ignore the possibility of eddy-induced Ekman pumping, a well known phenomenon where circulating ocean currents bring nutrients upwards or downwards within the eddies  [38, 39]. In this way we can isolate the plankton’s response to a single upwelling region, and study the effects of upwelling intermittency. Furthermore it allows us to simplify to horizontal advection only. In contrast to  [29], we analyse here a community that consists of two phytoplankton species competing for a limiting resource and grazed by zooplankton. The population model  [40] chosen displays excitability, which arises from the interplay of the fast dynamic timescale of phytoplankton growth (activator) with a slow development of zooplankton (the inhibitor). We chose two scenarios with different plankton communities: (I) where the community structure is shaped only by the availability of nutrients in the environment; (II) where both, the grazing pressure and the nutrient availability, trigger the bloom formation. First we show that these two systems exhibit very different spatio-temporal dominance patterns and display distinct and characteristic dynamical reactions to an upwelling event. Then we show that for both parameterizations even identical pulses of nutrient influx, trigger a diverse set of reactions in the plankton dynamics. The variety of possible responses can only be understood by analysing the interplay of different time scales that characterise the system as well as the interplay between the upwelling and mesoscale hydrodynamic structures present in the flow. The outcome is even more complex for the case of irregular pulses with a variety of strengths and durations. Our analysis shows that it is impossible to establish a relation between the HABs formation and upwelling events, by only looking at the respective time series of nutrients and plankton abundances without considering the mesoscale mixing by the ocean flow in the observation region.

The work is organised as follows: First, in Sec. 2 we briefly introduce the coupled hydrodynamic-biological model used. In Sec. 3.1, we examine how the position and initial time of nutrient parcels initialized at the upwelling region affect the residence time of nutrients in the observed area. In Sec. 3.2 we describe the spatio-temporal patterns of plankton for the scenario without upwelling and for the scenario subjected to a single upwelling pulse. Next, in Sec. 3.3 we describe how the response of the populations varies considering different initial times for an upwelling pulse. Furthermore, we analyse the chance of HAB formation for upwelling pulses of different duration and strengths. In Sec. 3.5 we extend our analysis to the study of the response to a series of irregular (intermittent) upwelling events. We conclude in Sec. 4.

2 Model

This section describes the modeling framework used in this work. The model consists of a two dimensional kinematic velocity field coupled to a biological model, see Fig. 1.

2.1 Hydrodynamic model

The hydrodynamic model is represented by a two dimensional kinematic velocity field $(u_{x},u_{y})$ that flows trough a predefined observation region passing by a circular obstacle (of radius $r$), located at $(x_{0},y_{0})=(0,0)$, mimicking an island. The flow velocity is such that it allows for the formation of vortices in the wake of the island. These vortices are released and carried away from the island along the observation region, from left to right in Fig.1 (a). Although this flow field can be obtained as a solution of the Navier-Stokes equation, we use an analytically generated field [41] that captures the main characteristics of this solution, but with smaller numerical efforts. The flow is characterized by a period T and in this approach a predefined stream function $\Phi$ is used to generate it: $u_{x}=\frac{\partial\Phi}{\partial y}$, $u_{y}=-\frac{\partial\Phi}{\partial x}$. This flow is known as an open chaotic flow in literature, for a detailed description of the modeling approach and the stream function see  [41, 29]. The model has been parameterized to represent one of the islands of the Canarian Archipelago, located in the Eastern Boundary Upwelling System off the African Coast, in agreement with [29], for details see Supplemental material.

2.2 Biological model

The biological model used, consists of a food web with three trophic levels NPPZ (Nutrients, two Phytoplankton species and Zooplankton) formulated in  [40] to describe the formation of harmful algal blooms (HABs). One of the phytoplankton species is considered to be toxic and the other one non-toxic, their concentrations are $P_{T}$ and $P_{N}$, respectively. They compete for a limiting nutrient resource, $N$, while being grazed by zooplankton, $Z$ (see Fig. 1 (b)). The inter- and intraspecific interactions are described by:

[TABLE]

The functional responses and parameterization are listed in the Supplemental material. As we will show below the system’s response to nutrient influx from upwelling results from a combination of bottom-up and top-down controls. To be able to analyse how these controls drive HAB formation we chose two different parameterizations for the population model: system (I), where the community structure and the dominant species results mainly from the availability of nutrients in the environment; system (II), where the grazing preference of zooplankton together with the nutrient availability both establish the resulting community structure. The parameters chosen for the two systems are very similar, with a few differences which emphasize different ecological processes. In system (I) the nutrient conversion rate $\theta_{i}$ and the respiration rates $r_{i}$ are different for each species ($\theta_{N}<\theta_{T},r_{N}>r_{T}$), so that the net growth rate of $P_{T}$ is always larger than the one of $P_{N}$ (see Fig. 2(a)). On the contrary, in system (II) these parameters are set to the same values ($\theta_{N}=\theta_{T},r_{N}=r_{T}$), for both species. Consequently, in system (II) the abundance of nutrients cannot drive a dominance change (see Fig. 2(b)). To especially test the role of grazing, we modify the parametrisation of zooplankton in the second system, in order to have a stronger influence of the grazer within the food chain. To that end, we boost the abundance of zooplankton by increasing its maximum grazing rate $\lambda$ and its growth efficiency on the non-toxic species $\alpha_{N}$. Finally we also add in this set up a strong preference of zooplankton to feed on non-toxic species ($\phi=0.05$), contrarily to system (I) where there is no preference ($\phi=0.5$). With this change, system (II) has a very strong top down control by design. This is demonstrated in Fig. 2 (c, d) that show an approximate net growth of zooplankton for the two systems. While for system (II) the net growth rate of zooplankton is always positive Fig. 2 (d), this does not hold for system (I). For very low abundances of phytoplankton, there is not enough food for zooplankton to survive. Therefore, for very low nutrient supply and subsequently very low phytoplankton abundance, zooplankton would go extinct.

As described previously the model takes into account the vertical influx of nutrients from the deep ocean into the mixed layer where all biological processes take place. The rate of this influx is given by $k$. This influx of nutrients may occur due to turbulent diffusion ($k=k_{d}$) or by vertical transport due to upwelling ($k=k_{\text{up}}$). The diffusive flux is $D_{v}\frac{\partial^{2}N}{\partial z^{2}}\sim D_{v}\frac{(N_{0}-N)}{h^{2}}$, where $N_{0}$ is the concentration of nutrients below the mixed layer and $h$ is an average extension of the gradient. By using the definition $k_{d}=\frac{D_{v}}{h^{2}}$ we can rewrite the relation as $k_{d}(N_{0}-N)$. For the ocean we find in the literature values of $D_{v}\sim 0.1-2.6$ m2 day*-1* [2, 42]. We use the known extension of the thermocline to estimate $h$ and therefore adopt values from $10$ to $25$ m. With these parameters we can evaluate $k_{d}$ in the range of $10^{-2}$ — $10^{-4}$ day*-1*. The nutrient transport due to upwelling, on the other hand, is defined as $u_{z}\frac{\partial N}{\partial z}\sim u_{z}\frac{\Delta N}{\Delta z}$. Therefore for the situation with upwelling we can define the thermocline exchange rate as $k_{\text{up}}(N_{0}-N)$ with $k_{\text{up}}=\frac{u_{z}}{h}$ + $k_{d}$. It is known that the vertical velocity $u_{z}$ may reach values as large as $\sim 40\;$ m day*-1* [42], however specifically for the region of the Canarian Archipelago we find estimations close to $10\;$ m day*-1* [43]. This gives us $k_{\text{up}}$ of 1 day*-1*. Therefore in this work we restrict ourselves to value for $k_{\text{up}}$ of the order of unity.

As already mentioned the nutrient influx is regulated by two parameters: the cross thermocline exchange rate $k$ and the nutrient concentration below the thermocline $N_{0}$. It is therefore compelling to outline here how the coupled effect of these parameters is reflected in the steady state community structure. These results are summarised in Fig. 3 (a) and (d) for system (I) and (II) respectively. The region of dominance of the toxic species is shown in red, and of the non-toxic in green. We also set-up two values for the cross thermocline exchange: $k_{d}$ and $k_{\text{up}}$ representing the conditions without and with upwelling respectively.

Our choice for system (I) corresponds to $k_{d}=0.015$ and $k_{\text{up}}=1.5$. These values will be used in all further simulations of system (I), unless stated otherwise. The dynamics of the system towards the steady state for diffusive exchange $k_{d}=0.015$ is shown in Fig. 3 (b) and leads to a steady state with a dominance of the non-toxic species, see Tab. 1. As explained previously the community structure in the system (I) directly reflects the low amount of nutrients of this scenario. Please note that for $k_{d}=0.015$ the presence of zooplankton allows for the coexistence of the two phytoplankton species. On the other hand, for a high nutrient supply $k=k_{\text{up}}=1.5$ the dynamics leads to the dominance of the toxic species and even the extinction of its competitor, see Fig. 3 (c).

The parameters chosen for the further analysis for system (II) are: $k_{d}=0.0045$ and $k_{\text{up}}=2.25$. Again we find, that with a large input of nutrients the toxic species dominates, see Tab. 1. However, the main mechanism how this dominance is achieved differs from previous case as explained previously. Please note that here we observe a significant amount of the total biomass concentrated at the higher trophic level, especially in the low nutrient limit. Another distinction of this set-up is the presence of both species of phytoplankton in the steady state for low as well as high nutrient influx Fig. 3 (e, f).

2.3 Coupled Model

The full biological-hydrodynamic model consists of the following reaction-advection-diffusion system of equations:

[TABLE]

where $C(x,y,t)$ represents the concentration of nutrients or plankton species in space and time, and $F_{c}$ are functions representing the biological interactions among these species, which are given by the population dynamical model Eq.(2.2). We consider a horizontal turbulent diffusion constant $D=10$ m${}^{2}s^{-1}$, that describes the advection by smaller scales in the flow field. Please note that Eq.(2.2) describes the dynamics as vertically averaged model only in the mixed layer while vertical transport is encapsulated in the biological model considering only the vertical exchange of nutrients and the sinking of phytoplankton (cf. subsection 2.2). The influx concentrations at the left boundary are setup as $20\%$ of the steady state values of Table 1. For numerical details please check Supplemental material. The code for the simulation reported here can be found in the Github repository: https://github.com/kseniaguseva/Upwelling.

3 Results

3.1 Hydrodynamic time scales

According to our aim to understand the conditions for HAB formation in the presence of an intermittent upwelling, it is important to analyze the interplay between hydrodynamic and biological time scales. We start the study of this nontrivial coupling by analysing the hydrodynamics that underlies all the biological processes in our system. To that extend we follow the motion of non-reacting fluid parcels passively transported by the flow field (tracers). Since we are interested in the impact of upwelling we compute the residence times of tracers starting in the upwelling region. Furthermore, we want to understand how the residence times of tracers depend on the initial time instant of their release. We measure it by releasing the tracers at a location ($x_{i},y_{i}$) and recording the time $\tau_{i}$ when they reach the right boundary at $x=8$.

We start by characterizing the possible trajectories that a tracer element can take depending on its release time $t_{i}$ and its release coordinates ($x_{i},y_{i}$), see typical trajectories and the respective $\tau_{i}$ in Fig. 4 (a). The main difference in $\tau_{i}$ arises from whether the tracer is captured by a vortex in the wake or not: The ones captured into a circular trajectory around the vortex core (black and gray trajectories in Fig. 4 (a)), spend at least two times longer in the observation area than the ones that are transported more or less straight by the main flow (red and blue trajectories of Fig. 4 (a)). Fig. 4 (b) summarizes our results on residence times of trajectories starting at ($x_{i}=0,y_{i},t_{i}$). The periodicity of the flow can be seen by the repeating patterns shown in Fig. 4 (b). The two finger-like structures in the residence times (blue points) correspond to tracers that are captured by the vortices. Another important characteristics of these patterns is their fractality (see Fig. 4(c)), which directly reflects the influence of the stable and unstable manifolds of the chaotic saddle present in the system  [41, 32]. Note that although we have chosen to fix $x_{i}$ at $x_{i}=0$, the results for other release positions within the upwelling region are very similar.

In summary, when we identify the tracers with nutrients released during upwelling, then the residence time of nutrients in the observation region changes with the position of the upwelling region, the location of release within that region, and the time instant when the nutrients are released. In other words, for how long nutrients, that have been released into the mixed layer during upwelling, are available for consumption by phytoplankton in the observation area depends crucially on the structure of the hydrodynamic flow at the time instant of upwelling. Next we will investigate the consequences that this effect has on the plankton dynamics.

3.2 Reaction to an upwelling pulse: spatio-temporal patterns

Before we start with the results of this section we shortly discuss the characteristics of the system in the absence of upwelling. In the absence of upwelling the spatial distribution of the biological species follow the uneven nutrient distribution in space. What is observed is an accumulation of nutrients in certain areas, the observed accumulation appears due to the small advective velocities around the island coupled with a constant nutrient influx from below the thermocline ($k_{d}$). Subsequently, this high nutrient concentration is captured by the vortices behind the island. This results into a bloom of non-toxic species in these regions in System (I), and a non-toxic bloom followed by the growth of zooplankton in System (II), see Supplemental material for details. We define the values of the spatial average over the observation region as $\left<C\right>$, where C stands for $N$, $P_{N}$, $P_{T}$ or $Z$. In addition, we will also use a distinct notation for the time average of $\left<C\right>$ for the case without upwelling, defining it as $\left<C\right>^{*}=\frac{1}{nT}\int_{0}^{nT}\left<C\right>dt$, the values for the two systems are shown in Table 2.

After having analysed the coupled hydrodynamic-biological model, let us characterize the HAB formed in the two systems in the presence of a simple upwelling pulse. We start by comparing the spatio-temporal patterns for the two biological systems for a case where the upwelling event triggers a HAB. We introduce a single upwelling event at $2.5$ T, at this instant the value $k_{d}$ at the upwelling area is exchanged to $k_{\text{up}}$. It is kept at this constant value for some time interval $\delta=0.5$, and then it is exchanged back to $k_{d}$.

In system (I) the spatio-temporal distribution is simple: the non-toxic species grows in the vortex cores undisturbed by the upwelling event, while the toxic species grows mainly by feeding on nutrients released by the upwelling. Fig. 5 illustrates these dominance patterns for three instances of time that follow an upwelling event. Here we have specifically chosen an upwelling event that triggers a strong dominance change. Please note the fast reaction of the toxic species to the nutrient influx.

By contrast, in system (II) both phytoplankton species readily grow in response to the nutrients. However there is a stronger response of the non-toxic species due to its lower half saturation constant, which allows it to reach high concentration and initiate the growth of zooplankton. However, the zooplankton development is a slow process and it only reaches significant concentrations when the non-toxic bloom is captured by a vortex. It is in this region where the toxic species, with extra nutrients and the presence of zooplankton, can successfully compete with non-toxic species. In fact, the high grazing pressure of zooplankton on the non-toxic species allows for the very localized dominance of the toxic specie, see Fig. 6. Note that when the bloom of the toxic species forms, the nutrients brought by the upwelling were already partially consumed.

We would like to emphasize that the spatio-temporal dominance patterns that appear in this system in the presence of upwelling in system (I) and (II) strongly differ. This difference can be explained by the fact, that the two spatio-temporal patterns result from distinct biological mechanisms. The behaviour in system (I) is solely determined by the bottom up control relying only on the supply of nutrients leading to a strong advantage of the toxic species in areas of high nutrient concentration. By contrast, in system (II) the top down control by the zooplankton is the dominant biological process shaping the spatio-temporal pattern. The toxic species can only dominate in areas, where its competitor is kept at low concentration due to the high grazing pressure. Furthermore, the two set-ups are characterised by different response times of the toxic species to the inflow of nutrients in the two systems.

3.3 The impact of initial time of the upwelling event

In Sec. 3.1 we have illustrated that fluid parcels released at different times, $t_{i}$, from the upwelling region can take very distinct paths trough the observation area. Some of these paths transport the fluid parcels directly away, describing a quick escape from the observation area, while others consists of spiral trajectories around vortex cores. These latter trajectories, in turn, are characterized by long residence times. In this section we will connect the advection with the plankton dynamics. Our objective is to answer how these different time scales affect the formation of HABs. Thus, we initialize upwelling pulses starting from different initial times $t^{\prime}$, with a predefined duration $\delta$ and strength $k_{\text{up}}$.

Now we analyse how these upwelling events impact the time series of the spatial averages of the plankton species of our biological model. While in the time series of system (I) only the toxic species exhibits a strong response to the upwelling events, in system (II) we observe, on the contrary, spikes in the growth of both phytoplankton populations and even in the abundance of zooplankton (Fig. 7 (a,d)). Despite these differences, we observe that in both systems the dynamics of the response of the plankton model to the upwelling event depends on its initial time $t^{\prime}$: in both systems we can have weak or strong responses, see Fig. 7 (a, d) and (b,e) respectively. This result is summarised in Fig. 7 (c,f) where the biomass of the toxic species is compared to the total biomass for an event with duration $\delta=0.5$. While for system (I) which is solely nutrient controlled we observe a dominance change for the average concentrations, this behaviour is absent for system (II), which has a strong top down control element. In system (II) we observe only local dominance change which never reaches a dominance of the toxic species in the spatial average. Please note the similarity of the diagrams of the two biological systems. The similarity of the response patterns for both systems (I) and (II) with respect to the timing of the response is entirely determined by the hydrodynamics.

3.4 Impact of the duration and strength of the upwelling

event

So far we have seen that the initial time of an upwelling event plays a crucial role for the mechanism of formation of HABs. In this section we extend our analysis to investigating the effect of the duration and the upwelling strength of randomly initialized upwelling events.

We compose the sequence of upwelling pulses in the following way: The upwelling events are initiated at particular time instants given by $4Tn+t^{\prime}_{n}$, where $n\in$ $\mathbb{N}^{*}$ and $t^{\prime}_{n}$ is chosen randomly for every $n$ from the interval $[0,T]$, see Fig. 8 (a). We establish that each sequence is characterized by upwelling events of duration $\delta$ and strength $k_{\text{up}}$. To quantify the effect of the upwelling on the growth of the phytoplankton species, the time series is divided into $n$ intervals: each one of them containing four periods and a single upwelling event. The time series of the average concentration $\left<C\right>$ for each one of these intervals is denoted $\left<C\right>_{n}$. Thus, the effect of each upwelling event on the population dynamics is reflected in the maximum, max$\left(\left<C\right>_{n}\right)$. Furthermore it is useful to systematically compare max$\left(\left<C\right>_{n}\right)$ to the average concentration in the absence of upwelling, $\left<C\right>^{*}$, we represent this deviation by $\Delta\left<C\right>_{n}=$ max$(\left<C\right>_{n})-\left<C\right>^{*}$ (see Fig. 8).

We start our analysis by fixing $k_{\text{up}}$. In the resulting time series of system (I), see Fig. 9 (a, b), the toxic species is the only species that shows a response to upwelling in its average values. On the contrary, in system (II) all the species show a bloom-like behaviour, Fig. 9 (c, d). It is clear that the average values shown for both of these systems, fail to completely describe the complexity of the spatio-temporal dynamics. Nevertheless, part of this complexity is revealed by the variability of different dynamic responses of the biological community to seemingly identical upwelling events, see Fig. 9. Comparing the different responses for the same system with the same duration, we notice that it depends crucially on the timing of the upwelling event, how strong the response is going to be. This revels clearly the importance of the structure of the flow field at the time instant of the upwelling. Additionally, our results reveal that this variability depends on the duration $\delta$, and this relation manifests itself in a similar way for both systems (Fig. 9). Our results reveal that longer upwelling events are associated with a vigorous growth of the toxic species. For this case the probability of HABs is large and we observe similar peaks in the concentration of the toxic species (large values of $\Delta\left<P_{T}\right>_{n}$). On the other hand, shorter $\delta$ values reveal a larger variety of possible outcomes. These results are summarized in the histograms of Fig. 10 (note the difference of the axes between the upper and the lower panels). The observed behaviour can be explained by taking into account that HAB formation depends on the temporal overlap between the upwelling event and the vortex formation in the wake of the island. Naturally for larger $\delta$ the probability of this overlap is higher and more nutrients are captured to incubate the growth of the toxic species. The strength of upwelling events has a complementary influence. For small values of $k_{\text{up}}$ the system needs longer upwelling events to release enough nutrients for toxic species bloom, see the Supplemental material for details.

At the end of this section we want to stress that from an analysis of the time series only, it is especially difficult to establish a causal relation between the upwelling event and the rate of increase of the toxic species. Although, an increase in the population of the toxic species always follows the upwelling in our model system, the level of increase in the population varies strongly, see Fig. 10 (d). This variety of the possible outcomes, however, can be easily explained by coupling the biological model with hydrodynamic mesoscale motion. Therefore by taking into account the interplay between the initial time of the upwelling event and the formation of vortices in the wake, it is possible to predict if the event will result in a HAB formation.

3.5 Intermittent upwelling events

In the previous sections we have seen that an upwelling pulse, even of the simplest possible profile, can result in a variety of possible outcomes for the plankton growth. The intricate interplay between plankton dynamics and the formation of vortices, or more general mesoscale hydrodynamic structures, results in time series showing responses of different strengths for identical upwelling events. Here in this section we analyse the response of our model to upwelling events that follow a time series that displays more complex patterns. The idea here is to mimic a more realistic situation, since upwelling is a wind driven phenomenon and hence, has an intermittent character. To generate this new time series of upwelling events we use a dynamical system which displays a special type of intermittent behaviour, known as “on-off” intermittency [44]. Two modes appear in this system: the “off” mode (situation without upwelling) where a very small value of an observable of the system sets up for long intervals of time; these intervals are interrupted by seemingly random bursts, characteristic to the “on” mode (upwelling events). Therefore the thermocline exchange rate at the upwelling region in the “off” mode is $k_{d}$ and in an “on” mode $k_{\text{up}}$, which here assumes a set of random values obtained from the dynamical system described in the Supplemental material.

Fig. 11 shows the response of our two biological systems to an identical sequence of intermittent upwelling pulses shown in light gray Fig. 11. Note that in this system there can be several short pulses of different strengths within a single period, furthermore the events are not isolated but come in small groups. Each group of pulses triggers a different outcome for the toxic population. Furthermore, as can be easily spotted in Fig.  11, there is a very strong variability between possible responses. Note, for instance the weak blooming behaviour of the toxic species around $t=2.5$ in contrast to the strong response at $t\sim 10$. For system (I) we find at $t\sim 2.5$ a rather long bloom with moderate amplitude, while at $t\sim 10$ the bloom exhibits a much higher amplitude. For system (II) we find a similar response, but now not only for the toxic but also for the non-toxic one and the zooplankton reflecting the importance of the grazing pressure in that case.

4 Conclusion

In this work, we have analysed how the competition, between two species for a shared limiting resource, can be affected by intermittent upwelling events providing an additional input of this resource. We have used a theoretical approach which couples the hydrodynamic flow field with a biological model by means of reaction-advection-diffusion equations. Notably, we have tracked the necessary environmental conditions that trigger a HAB. We were particularly interested on how the interplay of the hydrodynamic timescales as well as the mesoscale hydrodynamics structures, like vortices in the flow, coupled to intermittent upwelling pulses influence the spatio-temporal distribution of dominance patterns of different functional groups of phytoplankton. First we have characterized the HAB formation in two biological scenarios: the first scenario where the abundance of nutrients is the only factor responsible for the emergence of dominance patterns in the system; and the second one, where the dominance patterns arise from combination of competition for nutrients and grazing pressure from a higher trophic level. Both scenarios are characterized by distinct spatio-temporal inhomogeneous distributions of the phytoplankton groups, which appear as a result of an upwelling event. In the first scenario the toxic species develops along the whole nutrient plume, while in the second system a bloom is formed in a very localized region namely on a narrow ring around one of the vortices in the wake. The time of the bloom development also differs in these systems: in the first one the response of the toxic species is almost immediate, while in the second one the dominance change occurs while the vortex is advected away from the island. Despite the observed differences in theses two systems we demonstrate that, in both of them, the decisive factor triggering a bloom or not is the coupling of the upwelling event with the formation of mesoscale vortices. In this scenario the HAB formation results from the interplay of three timescales: (1) of the vortex formation at the island’s wake, (2) of the upwelling event and (3) of the biological growth. Our analysis shows that identical upwelling events that start in different instances of time may result in a variety of outcomes for the biological community depending on the properties of the flow at the moment of upwelling. The observed response depends on the time interval that nutrients released by upwelling spend in the observation region and consequently the quantity of nutrients captured by the vortices. Therefore the variability of these possible outcomes depends also on the duration and the strength of the upwelling events.

In summary we have observed that the HAB formation, independently of the biological set-up, is tightly associated with the transport dynamics of the flow field. From our analysis we conclude that without taking advection into account it appears to be not possible to establish the relationship between upwelling events and triggering a HAB. For this reason one cannot expect to find a functional dependence between upwelling events and plankton blooms in general, when only nutrients and plankton abundances are measured and no information about the flow field is available. Such measurements lacking the properties of the flow field will always be difficult to interpret and allow only conclusions when the flow field is simple and does not contain mesoscale hydrodynamics structures.

Acknowledgments

We are grateful to Rahel Vortmeyer-Kley for illuminating discussions. This work was supported by the Volkswagen Foundation (Grant No. 88459).

Supplemental Material

Hydrodynamic model

The analytically defined model describes the velocity field for an incompressible viscid fluid with a Reynolds number at which the solution of the Navier-Stokes equation is time periodic. The period of the flow is $T$. During this time, two vortices are created in the wake, with a phase shift of $T/2$, and move away from the island. The two vortices rotate in opposite directions and are characterized by a vortex strength $\omega$. One of them travels slightly above and the other slightly below the axis at $y_{0}$. Please note that the assumption of a two dimensional velocity field relies on the fact that the vertical velocities in the ocean are significantly smaller compared to the horizontal ones. Additional dynamical properties of the flow relevant to this work are reviewed in Sec. 3.1.

According to the situation in this geographical region the period of the flow $T$ is $32$ days. The parameters used for the flow field are shown in Table. 3. Also following  [29] we superimpose the Ekman flow $u_{E}$ in the $y$ direction, perpendicular to the main flow, for $x>1$, see Fig.1 (a) of the main text.

Biological model

The functional responses used and the parameters are listed in Table 4 and Table 5 respectively. Please note that these differential equations are based on some traditional NPZ models, such as of Steele & Henderson [45] and Edwards & Brindley [46]. An important characteristics of these models is that the nutrient uptake by phytoplankton $f_{N,T}(N)$ ($f_{N}(N)$ for the non-toxic species and $f_{T}(N)$ for the toxic species) is given by a Holling Type II functional response, while the grazing of zooplankton $h$ considers a Holling Type III functional response (see Table 4). Additional effects of interspecific and intraspecific competition are given by the function $g$, where $a/b$ is the maximum nutrient uptake rate of phytoplankton averaged over the depth of the mixed layer. The differences between the two groups of phytoplankton can be introduced through different parameters: their nutrient conversion rates $\theta_{N,T}$, half saturation constants $e_{N,T}$, respiration rates $r_{N,T}$, their feeding preference by zooplankton, $\phi$, and their quality as food for zooplankton expressed by the conversion rates $\alpha_{N,T}$. However, please note, that there is no direct influence of the toxic species on the mortality of zooplankton. Therefore this model is not restricted to HAB formation, but can be also used to describe the emergence of any phytoplankton bloom, in which the two different competing species are involved. The notation of toxic and non-toxic species simplifies the extension of our findings. The recycling by bacteria is considered indirectly with factors $\beta$ and $\gamma$ for conversion of the dead material back into nutrients.

Coupled Model

The full system of equations is solved using a semi-Lagrangian algorithm [31], the code can be found at https://github.com/kseniaguseva/Upwelling. We use a grid of $[500\times 300]$ points, the integration is carried out for step size $dt=0.01$, and the diffusion step $dt_{D}=\frac{dt}{10}$ which guaranties that $\frac{Ddt}{dxdy}<0.5$, so that the stability condition of the integrator is fulfilled.

Everywhere in the observation area except for the prescribed upwelling region located above the island, the cross thermocline exchange rate is set to $k_{d}$. In the upwelling region the value is exchanged between $k_{d}$ and $k_{\text{up}}$ in time, corresponding to intermittent upwelling events. The upwelling region, if not stated otherwise, spans the region: $x\in[-1,1]$ and $y\in[2,2.5]$, see Fig.1 of the main text.

All the modeled species enter the system from the left at $x=-2$ with the same concentration at all $y$ values and are advected across the observation area. We assume that they arrive from the open ocean, an environment poor in nutrients and plankton. Therefore we use as the influx $20\%$ of the steady state concentration value reached by each given species for a cross thermocline exchange rate $k_{d}$. Please note that all species of phytoplankton are present in the system in the influx. As we show in Sec.2.2 of the main text the non-toxic species dominates in the influx conditions, since those are low in nutrients, for both scenarios that we analyse.

Biological model with hydrodynamics, in the absence of upwelling

Hare we present the dynamics of the coupled biological-hydrodynamic model without upwelling to allow a comparison with the results in the main text. First, we can mention that system (I) has very similar properties as the population dynamics analysed in [32]. But in contrast to [32] we are here more interested in the dynamics of plankton in the whole area and not only in the development of a plankton bloom related to specific regions such as vortices. When the upwelling is not present the non-toxic species is found in higher concentrations in a region around the island and inside every vortex formed in the wake. This spatial distribution of the non-toxic species reflects the nutrient accumulation in those regions, which results from the vertical exchange of nutrients across the thermocline while the horizontal advection is slow. Subsequently, this high nutrient concentration is captured by the vortices behind the island, where it creates good conditions for the growth of the non-toxic species. Nevertheless this accumulation of nutrients does not reach the threshold to trigger a HAB. Notably, system (I) is characterized by a periodic timeseries for the spatial average of the non-toxic species following the vortex formation and advection: $\left<P_{N}\right>$ has a period $T/2$. For system (I) the averages correspond to: $\left<P_{N}\right>^{*}=0.042$ gC m*-3*, $\left<P_{T}\right>^{*}=0.019$ gC m*-3* and $\left<Z\right>^{*}=0.005$ gC m*-3*, $\left<N\right>^{*}=0.025$ gC m*-3*.

We continue with the characterization of the spatio-temporal patterns formed in system (II). For the scenario without upwelling the non-toxic species, as in the previous case, develops around the island and inside the vortices. As in system (I), the growth of phytoplankton also follows the large concentration of nutrients. The large presence of non-toxic species creates now good conditions for the development of zooplankton, that also grows in the same location but at a slower rate. As a result, this system displays more complex spatio-temporal patterns, where as the vortices are advected by the flow field the concentration of zooplankton increases as well. On the other hand due to the large phytoplankton growth these vortices become depleted in nutrients. Therefore, without upwelling, conditions that would give an advantage to the development of the toxic species are never met. In summary, in system (II) the non-toxic species is the dominant species everywhere in space with the average concentration $\left<P_{N}\right>^{*}=0.0084$ gC m*-3*, which is in agreement with low nutrient concentration in the system $\left<N\right>^{*}=0.0286$. The average concentration of the other plankton species are: $\left<P_{T}\right>^{*}=0.0010$ gC m*-3* and $\left<Z\right>^{*}=0.0054$ gC m*-3*. Furthermore, the average concentration of all species is about one order of magnitude smaller than for system (I). This is a direct consequence of the used values of the thermocline exchange rate.

Nutrient release from the upwelling region

Finally, we investigate the last parameter that controls the release of nutrients from the upwelling region — the cross thermocline exchange rate $k_{\text{up}}$. To get some insight of the effect of $k_{\text{up}}$ on the growth of toxic species we compare the $\Delta\left<P_{T}\right>_{n}$ to the undisturbed concentration value $\left<P_{T}\right>^{*}$. More precisely we calculate $\frac{1}{n}\sum_{n}\frac{\Delta\left<P_{T}\right>_{n}}{\left<P_{T}\right>^{*}}$, for $n=150$ upwelling events for a range of $k_{\text{up}}$ and $\delta$ values. As expected small values of either $k_{\text{up}}$ or $\delta$ do not trigger the growth of the toxic species. By contrast, for high values of these parameters, the toxic species exhibits high growth in both systems. Fig. 12 shows the impact of these two parameters simultaneously. Please note that $k_{\text{up}}$ may limit the HAB formation given a fixed $\delta$. The mechanism behind that is simple: the smaller $k_{\text{up}}$, the more time is needed for the released nutrients to reach a significant level for the toxic species development. This means that for a small value of $k_{\text{up}}$ the initial times of upwelling events are even more constrained, i.e. they have to start at the beginning of vortex formation to allow for sufficient input of nutrients. In Fig. 12 we also show how the dominance change of species depends on parameters $\delta$ and $k_{\text{up}}$. To this end we compute the probability, $P_{\text{HAB}}$, that the toxic species out-competes the non-toxic one as a result of an upwelling event, $P_{\text{HAB}}=P(\text{max}(\left<P_{T}\right>)_{n}>\text{max}(\left<P_{N}\right>_{n}))$. This probability again is computed using the time series with $200$ upwelling events. In Fig. 12 (a) we delineate the region $A$, where $P_{\text{HAB}}>0.98$. In this region most of the upwelling events lead to HABs. The region $C$ is restricted to parameters where $P_{\text{HAB}}<10\%$, while the intermediate region $B$ spreads over the parameter range where many different outcomes are possible. Note that the sizes of these regions depend on the parameters of the biological model. Since in the system (I) mainly the toxic species responds to the upwelling event, the largest part of the parameter spaces is covered by the region A. By contrast, as we have explained before, the dominance change in system (II) is restricted to some spatial regions and therefore, region $C$ spans over the whole parameter space for this scenario, see Fig. 12 (b).

Appendix A Intermittency

To model the time series of intermittent upwelling events we use the absolute value of the $x_{1}$ variable from the system of the following differential equations:

[TABLE]

where $\mu=1.815$, $\nu_{1}=1.$, $\nu_{2}=1.815$, $\nu_{3}=0.44$, $\nu_{4}=2.86$, $\nu_{5}=2.86$. For these parameters the system displays on-off intermittency, for details see [44].

To lead to upwelling events of the adequate duration we rescale the time in Eqs.(3) using $t^{\prime}=t/T$, where T is the period of the flow field. The “off” states of the time series from Eqs.(3) are characterized by $|x_{1}|\sim 0$, and the “on” state by max$(|x_{1}|)\sim 1$. Therefore, we transform this time series to:

[TABLE]

With this transformation the new time series has the “ off“ state $\sim k_{d}$ and an “on” state which can be at most $k_{\text{up}}$. The time series $k(t^{\prime})$ is used to define the cross thermocline exchange rate at the upwelling region.