Impact of critical eddy diffusivity on seasonal bloom dynamics of Phytoplankton in a global set of aquatic environments

The intensity of eddy diffusivity and the spatial average of water velocity at the depths of the water column in oceans and lakes play a fundamental role in phytoplankton production and phytoplankton and zooplankton biomass, and community composition. The critical depth and intensity of turbulent mixing within the water column profoundly affect phytoplankton biomass, which depends on the sinking characteristic of planktonic algal species. We propose an Nutrient-Phytoplankton-Zooplankton (NPZ) model in 3D space with light and nutrient-limited growth in a micro-scale ecological study. To incorporate micro-scale observation of phytoplankton intermittency in bloom mechanism in stationary as well as oceanic turbulent flows, a moment closure method has been applied in this study. Experimental observations imply that an increase in turbulence is sometimes ecologically advantageous for non-motile planktonic algae. How do we ensure whether there will be a bloom cycle or whether there can be any bloom at all when the existing phytoplankton group is buoyant, heavier, motile, or non-motile? To address these questions, we have explored the effects of critical depth, the intensity of eddy diffusivity, spatial average of water velocity, on the concentration as well as horizontal and vertical distribution of phytoplankton and zooplankton biomass using a mathematical model and moment closure technique. We quantify a critical threshold value of eddy diffusivity and the spatial average of water velocity and observe the corresponding changes in the phytoplankton bloom dynamics. Our results highlight the importance of eddy diffusivity and the spatial average of water velocity on seasonal bloom dynamics and also mimic different real-life bloom scenarios in Mikawa Bay (Japan), Tokyo Bay (Japan), Arakawa River (Japan), the Baltic Sea, the North Atlantic Ocean, Gulf Alaska, the North Arabian Sea, the Cantabrian Sea, Lake Nieuwe Meer (Netherlands) and several shallower lakes.

The impact of the physical mixing process on the interaction of phytoplankton species is immense [1][2][3] .The occurrence of phytoplankton bloom gets triggered by a rapid increment in the spatial distribution of phytoplankton density.Fundamentally speaking, regional and seasonal changes influence environmental factors like solar irradiance and water salinity, which cause irregularities in phytoplankton growth.The nature of dominating or existing plankton classes changes with regional and seasonal changes.Different phytoplankton classes possess different traits, most notably size, shape, specific rate of growth, and behavior, such as the nature of floating or sinking, which all together determine their dynamical behavior in an aquatic ecosystem.
Eddy diffusion exerts an intriguing impact on the spatial distribution of plankton species along horizontal and vertical directions for any turbulent aquatic environment 4 .Changes in turbulent mixing cause a significant shift in the existing phytoplankton community.Increase in the oceanic mixed layer depth negatively affects the leading light available to planktonic algae and the sedimentation loss rates of sinking phytoplankton, resulting in an overall decrease in phytoplankton growth with increasing mixing depth 5 .While buoyant species mostly float upwards during weak mixing, intense mixing causes species to dive into a deeper zone.However, significant change in the intensity of vertical mixing triggers the transportation of nutrients and other essential ingredients required for photosynthesis from the surface layer to a higher depth.Nutrient enrichment positively affects nutrient availability and phytoplankton growth, resulting in persistence of sinking species into deeper zones, which initiates the chances of bloom development at higher depths.The nature of bloom gets dominated by the nature

Model definition: conventional nutrient-phytoplankton-zooplankton (NPZ) model
According to several observations regarding phytoplankton bloom occurrences in a global set of aquatic environments (Table 1), we initially propose an NPZ model in 3D space considering the contribution of turbulence mixing and the nature of phytoplankton species, where we have later applied the moment closure method for better understanding of bloom mechanism.Starting with non-averaged small-scale equations and ignoring molecular and turbulent diffusion, we first consider the space-time evolution equation for the NPZ ecosystem model, given by In this model, u, v are the horizontal water velocities, w is the vertical water velocity, k n , k p , k z are the hori- zontal swimming speed of nutrient, phytoplankton and zooplankton, w n , w p , w z are the vertical sinking velocity (1) The term (m P + r P + s + k 2 ) describes the overall phytoplankton loss due to normal mortality ( m P ), respira- tion ( r P ), sinking (s) and cross thermocline exchange ( k 2 ).The term g z PZ is the loss of phytoplankton due to zooplankton grazing, where g z is the grazing rate.Suppose a z is the food conversion efficiency of zooplankton.In that case, zooplankton growth includes a z g z PZ caused by assimilated grazing impact, whereas zooplankton loss includes excretion ( τ z Z ) and mortality due to competition ( m z Z 2 ).As the proposed system is conserved, there is no external source for nutrient N. The growth of phytoplankton causes nutrient loss, and the system regains nutrients from the loss of phytoplankton and zooplankton.

NPZ closure model
During the studies of the distribution of microscale phytoplankton using high-resolution profiling fluorometers, it has been observed that the local fluorescence values are highly fluctuating in space 7,18 .Since the density of zooplankton is interrelated to the density of phytoplankton and phytoplankton growth is interconnected to the density of nutrients available in the closed environment, all three variables N, P, Z are expected to be fluctuating in space in the closed aquatic ecosystem.However, the non-closure model as described by Eqs. ( 1)-(3) does S 1 (N, P, Z) = − µ p P + (m P + r P + s + k 2 ) + (1 − a z )g z PZ + (τ z + m z Z)Z, S 2 (N, P, Z) =µ p P − (m P + r P + s + k 2 ) P − g z P Z, S 3 (N, P, Z) =a z g z PZ − (τ z + m z Z) Z,  not consider this spatial variability.To derive the NPZ-closure model from the conventional NPZ model (nonclosure), the model variables are considered to be a function of both time (t) and space (r), namely, where N 0 , P 0 , Z 0 are spatial mean values of the nutrient, phytoplankton, and zooplankton communities respec- tively, and N ′ , P ′ , Z ′ are their respective fluctuating components corresponding to each mean value.The hori- zontal and vertical sampling for microscale phytoplankton distribution has the same statistics at the centimeter scale 19 and also at the millimeter scale (except for extreme values).Hence, the statistics of the fluctuating components is independent of the direction of sampling (isotropic).Therefore, the spatial average of each fluctuating component is zero at any particular time, that is, < N ′ (r) > = 0, < P ′ (r) > = 0, < Z ′ (r) > = 0, while its temporal average cannot be zero, which implies < N(t) >= N 0 (t) , < P(t) >= P 0 (t) and < Z(t) >= Z 0 (t) .Also, water velocities u, v, w are fluctuating components in space and time, whereas k n , w n , k p , w p , k z and w z are non fluc- tuating quantities, that is,

Modelling framework
Substituting (4-5) in (1), and applying the Reynold's averaging method in space, we obtain, since u , v , w , k n , w n are independent of x, y, z.
Three terms on the left hand side of the above equation namely, �u ′ N ′ �, �v ′ N ′ �, �w ′ N ′ � , are covariance terms between velocity fluctuations and fluctuations of nutrient density, and therefore denote turbulent transports.These are classically parameterized by means of a down-gradient approach as with horizontal and vertical eddy diffusion coefficient or turbulent diffusion coefficient or eddy diffusivity k h and k v respectively.Therefore, substituting these values of �u ′ N ′ �, �v ′ N ′ �, �w ′ N ′ � , we obtain Considering k h , k v to be constant parameters, we obtain, To construct a closure for these quantities, we follow the principle of second order turbulence closure modelling 20 and first derive transport equations for N ′ , P ′ and Z ′ by subtracting the equation for N 0 , P 0 , Z 0 from the equation for N and accordingly for P and Z. Once these equations are derived, we note that The consequent application of this derivation would result in vertical turbulent flux divergence formulations for the second order moments �N ′2 �, �P ′2 �, �Z ′2 �, �N ′ P ′ �, �P ′ Z ′ �, �N ′ Z ′ � , which are not in gradient form.Therefore, we ignore vertical turbulent fluxes in the derivation and later parameterize the vertical flux in the down-gradient form with the same eddy diffusivity as for the other tracers 21 , and obtain (4) N(r, t) = N 0 (r, t) + N ′ (r, t), P(r, t) = P 0 (r, t) + P ′ (r, t) and Z(r, t) = Z 0 (r, t) + Z ′ (r, t), www.nature.com/scientificreports/Using (6), and assuming that the variances and the covariances, �P ′2 � , �N ′2 � , �Z ′2 � , �N ′ P ′ � , �P ′ Z ′ � , �N ′ Z ′ � , fol- low the same diffusion law as N 0 , P 0 , Z 0 and assuming the covariances �N ′ P ′ �, �P ′ Z ′ �, �N ′ Z ′ � have horizontal and vertical settling velocity k np , k pz , k nz and w np , w pz , w nz respectively, and re-introducing vertical flux divergences for the second-moment equations, we finally arrive at the following set of equations: We have assumed that the random variables N, P and Z follow a joint lognormal probability distribution whose observed values are N, P, Z (densities of nutrient, phytoplankton and zooplankton respectively).The lognormal distribution can be fitted well to empirical data and has been widely used in continuous model 5 .We have ignored the third and higher order fluctuating terms to obtain simple closure.
Equations (10-12) represent the time evolution of mean terms, Eqs.(13-15) represent time evolution of variance terms and Eqs.(16-18) give the time evolution of covariance terms.Now the corresponding bio-physical model is ( 9) Vol:.(1234567890) www.nature.com/scientificreports/Summing up first three equations of the above system provides since 3 i=1 �S i � = � 3 i=1 S i � = 0 .Adding next three and twice of last three equations of ODE system, we obtain where B is the variance of the sum of all fluctuating components.Therefore, N 0 + P 0 + Z 0 and With appropriate scaling, the nine Eqs. (10-18) can be reduced to seven equations of dimensionless parameters and variables as where V = (p 0 , n 0 , �p ′2 �, �p ′2 �, �z ′2 �, �n ′2 �, �p ′ z ′ �, �n ′ p ′ �) and f 1 , f 2 , .., f 7 are functional forms of V provided in sup- plementary information (SI).All dimensionless parameters and variables are defined in Table 2. Clearly,

Model analysis
Let, V * = (p 0 * , z 0 * , �p ′2 � * , �z ′2 � * , �n ′2 � * , �p ′ z ′ � * , �n ′ p ′ � * ) be a spatially homogeneous steady state of dimensionless reaction-advection-diffusion closure system.Then, , where A is total biomass of the system Vol We define, to be a small inhomogeneous perturbations about the steady state V * .Provided these perturbations are suffi- ciently small, it is again possible to linearize the dimensionless reaction-advection-diffusion closure system about the homogeneous steady state V * .Now the Jacobian of the linearized system about the steady state is  29 w p m day −1 0.012-6 30 - www.nature.com/scientificreports/Hence, the linearized equations of the closure system can be written as We now substitute, Q i = α i e − t e −σ 1 x e −σ 2 y e −σ 3 z , , σ 1 , σ 2 , σ 3 > 0 in the above system and obtain 7 equations in α i , which can be written in matrix form as where 0 is zero vector of order 7, α = (α 1 , ..., α 7 ) T is a column vector and S is matrix of order × 7 , where all non-diagonal elements are same for the matrices S and the Jacobian (J) of the linearized closure system about the point V * .The diagonal elements s ii are as follows: S × α = 0, We now obtain the characteristic equation of the linearized closure system as Det(S) = 0 for some suitable choices of σ 1 , σ 2 , σ 3 and α = 0.

Analysis of critical threshold value
We now consider a particular set of parameter values (to be estimated) representing a particular circumstantial scenario and a suitable set for amplitudes σ 1 , σ 2 , σ 3 .For aquatic ecosystem, k p , w p are swimming and sinking velocities of phytoplankton community, which only depend on the nature of existing or dominating phytoplankton species, hence they do not vary with seasonal or environmental factors.Therefore, among all advection and diffusion factors, k h , k v , u , v , w are the control parameters.Hence, depending on the estimated parameter set, we consider the dominating control parameter to obtain its corresponding critical value.
Let us assume that, for a considered set, k v , w are the dominating control parameters of the system.Keeping k v as a variable and putting other values, we determine S. Thus, Det(S) = 0 provides the characteristic equation of degree 7, where each coefficient is a polynomial of k h .Here, Our aim is to construct a Routh Hurwitz matrix from this polynomial.Degree of the characteristic polynomial is 7, therefore dimension of the Routh-Hurwitz matrix will be 7 × 4 .Elements of first two rows are coefficient of the characteristic polynomial, other elements are determined by cross multiplying the coefficients of previous rows.
We now proceed with Routh Hurwitz stability criteria to find a domain of k v , which will cause the spatially uniform steady state to be stable.The elements of Routh-Hurwitz matrix are polynomial functions of k v and nature of stability of spatially uniform steady state along with the system depends on the nature of eigenvalues, which gets decided by the sign of elements of first column of Routh Hurwitz matrix.This provides a lower bound on the domain of k v , say, k v

Boundary and initial conditions
To conserve the property of closure model, zero flux boundary condition has been introduced to the system.While non-dimensionalizing the system we have considered h is the length scale along x, y directions and l is the depth scale along z direction, in particular, we have considered a particular cubic part of an aquatic system having volume h 2 l meter 3 .

According to zero-flux boundary condition, rate of change of model variables along vertical and horizontal directions remains zero at boundary points implying
The remaining variables z 0 , �p ′2 �, �z ′2 �, �n ′2 �, �n ′ p ′ �, �p ′ z ′ � satisfy the same conditions.The initial conditions are In the formulated model variables, all mean and variance terms are always non-negative, whereas covariance terms can take any sign (depending on the relation between the corresponding model variables) and, corresponding dimensionless variables also follow the same.These dimensionless terms satisfy the condition, Hence, distribution of these variables are bounded by 1.On the other hand, phytoplankton mean decreases with depth, depending on that, mean of herbivorous zooplankton varies, both are bounded by constant 1 ( p 0 + n 0 + z 0 = 1 ) for dimensionless system and all these dimensionless quantities p 0 , n 0 , z 0 are strictly non-negative.Therefore, the functions φ 1 , ..., φ 7 are chosen in such a way so that they satisfy these criteria.
Besides, in this model it is assumed that only phytoplankton community has swimming and sinking velocities, hence k p , w p > 0 , whereas k n = 0, w n = 0, k z = 0, w z = 0 .Also, we have assumed that k np = 0, w np = 0 , k pz = 0, w pz = 0 (since we lack any observational or theoretical basis for modelling the effect of swimming and sinking on covariance).

Critical depth ( D c )
The phytoplankton species existing in higher depth should have access to essential components like solar irradiance for the continuation of photosynthesis, and the availability of such ingredients depends on the value of the depth (say, H c ) at which sinking phytoplankton species get stuck because of the resilience of water body.For any water column, there exist a layer (say, D c ) beyond which total loss of phytoplankton biomass neutralizes the total growth of phytoplankton, so bloom cannot occur underneath this depth D c .This depth is called critical depth.Depending on the value of this depth, the possibility of bloom occurrence in higher depth is decided.Also, depending on the fact whether H c is less than D c or not, the possibilities of the occurrence of bloom in higher depth within the water column vary.Critical depth ( D c ) is calculated using a reformation of the Sverdrup equation 22 , namely, where E 0 is the surface photosynthetically active radiation (PAR) integrated over 24h expressed in mol photons m −2 day −1 and κ is light attenuation coefficient.

Results
To capture the role of physical mixing process and spatial average of water velocity on bloom dynamics in global scale, our numerical investigation (the assigned parameter values are obtained from Table2) is sectioned into two parts, (i) System 1: spatial average of water velocity is zero ( �u� = �v� = �w� = 0 ), but turbulent diffusion is non- zero ( k h = 0, k v = 0 ), which is possible when mixing of comparatively highly dense or less dense external fluid or temperature difference causing a difference in water density which gives rise to diffusivity without influencing water velocity, like the lake, pond and almost calm marine ecosystem; (ii) System 2: spatial average of water velocity along with turbulent diffusion is non-zero ( �u� � = 0, �v� � = 0, �w� � = 0, k h � = 0, k v � = 0 ), which stands for a complete turbulent medium like rest- less marine environment and flowing river.

Surface bloom (SB)
SB1: When phytoplankton community is buoyant in nature ( k p > w p ) and spatial average of water velocity is zero ( �u� = �v� = �w� = 0) Figure 1a represents the vertical distribution of phytoplankton biomass of buoyant community ( k p > w p ) on the 60th day of the warmer season ( I 0 = 320 W/m 2 ) through a water column having depth L z = 60 meters in where γ = σ 3 2 and C i , i = 1, ..., 7 are non-negative constants.an aquatic medium of system 1.In contrast, vertical turbulent diffusion ( k v ) is low as high surface irradiance I 0 = 320 W/m 2 causes high thermal stratification of surface water resulting in reduced vertical eddy diffusivity.However, it dominates horizontal eddy diffusivity ( k v > k h ) in the absence of surface wind.Under such cir- cumstances, the surface bloom of the buoyant community has been observed on the 60th day of warmer season while existing low vertical turbulence is highly dominated by its corresponding critical value, that is, k v < k v c . Figure 1b represents the horizontal distribution of phytoplankton during that period, where we observe several yellow to dark green patches floating on surface water, which indicates intermittent patches of phytoplankton bloom covering the surface water.Here, the numerical results show that high buoyancy, in the presence of low grazing influence, exerts a positive impact on phytoplankton productivity and growth rate V m , by driving domi- nating phytoplankton community to light and nutrient-rich zone (Fig. 1a,b).SB2: When phytoplankton community has higher sinking tendency ( k p < w p ) and spatial average of water velocity is zero Figure 1c represents the vertical distribution of phytoplankton biomass of heavier community on the 115th day of the warmer season ( I 0 = 220 W/m 2 ) through a water column having depth L z = 60 meters in an aquatic medium of system 1.In contrast, vertical eddy diffusivity ( k v ) is low, which still dominates horizontal eddy diffusivity ( k v > k h ) during low wind mixing period in the presence of higher grazing pressure.Under such circumstances, a surface bloom of the existing phytoplankton community has been observed on the 115th day while critical vertical eddy diffusivity dominates existing vertical eddy diffusivity ( k v < k v c ) and zooplankton community is abundant on the surface zone (Fig. 2a).The figures (Figs. 1c, 2b) indicate that higher phytoplankton biomass provides ample food to grazers (zooplankton community), hence this aquatic food chain relation generates the possibility of bloom development of both of the species on the surface layer.SB3: When phytoplankton community is buoyant in nature ( k p > w p ) and spatial average of water velocity is non-zero ( �u� � = 0, �v� � = 0, �w� � = 0) Figure 1d stands for the occurrence of surface bloom of buoyant species in an aquatic medium of system 2, in the presence of weak vertical eddy diffusivity ( k h > k v ) in an intense wind mixing period during the warmer season ( and grazing pressure does not hinder phytoplankton growth on surface water (Fig. 2a).The figure shows the vertical distribution of phytoplankton biomass ( p 0 ) on , (g) Vertical distribution of p 0 of heavier community of system 2 ( k v > k v c , w > w c ) (h) Vertical distribution of p 0 of buoyant community of system 2 ( k v > k v c , w > w c ).All the figures are in the presence of low grazing influence.The parameters values are mentioned in SI. the 39th day of the warmer season (late spring/ summer/ early autumn) through a water column having a depth of 60 meters.Figure 1d represents the vertical distribution of phytoplankton biomass ( p 0 ) of the same com- munity in a turbulent medium through a water column having a depth of 60 meters on the 39th day of summer season ( caused by the eventual stopping of the wind event.Therefore, in the presence of low grazing influence, intense wind mixing period on the horizontal layer gives rise to the fact , which initiates the process of surface bloom formation of phytoplankton community (Fig. 1e).In contrast, the eventual stopping of wind event under the same grazing pressure causes k h < k h c , u < u c , v < v c and hinders surface bloom (Fig. 1d) and depending on the nature critical vertical eddy diffusivity k v c and the spatial average of water velocity w , nature of bloom formation or phytoplankton biomass accumulation at different layers varies accordingly (Fig. 1d).SB4: When phytoplankton community is heavier ( k p < w p ) in nature and spatial average of water velocity is non-zero Figure 1f represents vertical distribution of phytoplankton biomass ( p 0 ) of heavier species with higher sinking tendency in the presence of low grazing influence and weak vertical eddy diffusivity ( , w < w c on 81th day of warmer season.Eventual stopping of the wind event during a thermally stratified period gives rise to the fact , which activates transportation of sinking biomass to upper layers in the absence of wind event and generates the possibility of bloom formation on the upper ocean in the presence of low grazing pressure.

Bloom in higher depth (BHD)
BHD1: Phytoplankton community has higher sinking tendency ( w p > k p ) Figure 1g represents the vertical distribution of phytoplankton biomass on the 50th day of the cold season ( I 0 = 160 W/m 2 ) through a water column having depth L z = 60 meters in an aquatic medium of system 2, where the phytoplankton community has higher sinking tendency ( w p > k p ) in the presence of high vertical eddy dif- fusivity satisfying k v > k h , since low surface irradiance I 0 = 160 causes less thermal stratification of surface water which triggers kinetic energy of water particles, hence vertical turbulent diffusion k v increases.It dominates the horizontal eddy diffusivity k h in the absence of a wind event.The intra-specific competition among herbivorous ), (c) Vertical distribution of p 0 on 60th day of cold season ( k v < k v c , < w >< w c ) of system 2, p 0 has higher sinking tendency ( w p > k p ), (d) Vertical distribution of p 0 on 48th day of cold season of system 1 ( k v << k v c ), p 0 has higher sinking tendency, (e) Horizontal distribution of p 0 at the depth H c = 60 meters on 65th day of cold season ( k v > k v c , < w >> w c ) of system 2, p 0 has higher sinking tendency, (f) Relation between w and k v c using estimated data set, (g) Vertical distribution of p 0 of buoyant community ( k p > w p ) on 35th day of cold season of system 2 ( k v < k v c , w < w c ), (h) Vertical distribution of n 0 of system 2 on 32nd day of cold season ( k v < k v c , w < w c ), p 0 is buoyant in nature ( k p > w p ).The parameters values are mentioned in SI. www.nature.com/scientificreports/zooplankton classes for the same food sources remains high, and the grazing and assimilated grazing coefficient remains low.This is the case where a lack of thermal stratification triggers vertical eddy diffusivity to initiate the possibility of phytoplankton bloom formation at higher depths in the presence of low grazing pressure.
BHD2: When phytoplankton community is buoyant in nature ( w p < k p ) Figure 1h represents the vertical distribution of phytoplankton biomass of buoyant community on the 70th day of the cold season ( I 0 = 160 ) through a water column having depth L z = 60 meter in the presence of high vertical eddy diffusivity ( k v ) satisfying k v > k h (which happens when no wind event affects surface waves of turbulent flow) in a turbulent medium ( �u� � = 0, �v� � = 0, �w� � = 0, where w > u , v ).According to this figure, when critical vertical eddy diffusivity is crossed ( k v c < k v ) and vertical water advection dominate its critical value ( w > w c ) in the presence of low grazing influence ( g z = 0.4 ), species sink to higher depth and gets accumulated at the bottom layer ( L z = 60 < D c , SI Table 2) and eventually bloom occurs over there.
Note: Relation between controlling factors and depths are provided in SI Tables 1, 2 of system 1 and system 2 respectively.

Discussion
Phytoplankton experience turbulence as an instantaneous, linearly varying fluid velocity across the cell body 32,33 , and biomass distribution 34 .Nonmotile phytoplankton such as some diatom and cyanobacteria species orient in hydrodynamically favorable directions to displace over large distances 35 by either positively or negatively varying their buoyancy through different mechanisms such as gas vesicles, lipid storage, the exchange of heavy ions with their surroundings or by forming chains that modify their sinking rates [36][37][38] .Even though it is energetically demanding, motility is advantageous in a world where resources (nutrients) are heterogeneously distributed 39 .The majority of the phytoplankton classes are flagellated and, therefore, motile.Marine pelagic plankton such as ciliates, flagellates, and copepods exhibit a wide variety of motility patterns such as straight line and helical swimming, darting motion, etc.Cells move to new locations to seek more favorable environments or to interact with other organisms.Along with other variety of cues, including their encounter with other individuals (prey, a predator, or a mate), or environmental conditions such as light 40 , temperature, etc., these complex motility patterns are behavioral responses of plankton, which is controlled by turbulence 41 , which vary with regional and seasonal changes.Water turbulence significantly affects the distribution of phytoplankton 42 , and this leads to the concept that the formation of phytoplankton patches depends on the balance between phytoplankton growth rate and nature of eddy diffusivity, which controls the diffusion of the water 43 .
For any aquatic medium, the occurrence of bloom depends on how eddy diffusivity drives the distribution and dynamics of existing communities and how the existing community successfully dominates overcoming all hindrances.In order to investigate the reasons behind such occurrences with seasonal variation, we have focused on five components in our numerical investigation, namely, (i) nature of dominating phytoplankton community, (ii) horizontal and vertical eddy diffusivity, (iii) spatial average of water velocity, (iv) critical depth ( D c ), and (v) grazing impact.
We now consider the case of occurrence of bloom at a higher depth H c for predefined systems considering two scenarios, (i) when vertical eddy diffusivity ( k v ) is the only dominating component in system 1; (ii) when the vertical component of the spatial average of water velocity dominates over the corresponding horizontal component ( w > u , v ) and vertical eddy diffusivity dominates horizontal eddy diffusivity ( k v > k h ) in system 2.
We now address seven major research questions, which are the focus of our work.
The first question which we like to address is why phytoplankton bloom in higher depth is so common in nature during the winter season, and when can bloom also be developed in the upper ocean or surface water in the same season?
During the cold season, due to less thermal stratification of surface water for both predefined systems, the stability of the water column gets disturbed, causing an increment in vertical turbulent diffusion in the less stratified layer.Besides, externally induced artificial mixing causes an imbalance in the vertical velocity distribution of stationary flow and generates a turbulent environment, which triggers vertical eddy diffusivity k v .As k v gets affected, vertical mixing gets affected as well; as a result, the process of stabilization gets started to restore the system's stability through the transportation of biomass carried within the water column.While this stabilization process continues, phytoplankton biomass gets shifted from a highly concentrated zone to a less dense zone.Under such circumstances, considering the grazing impact to be low, several scenarios may arise regarding the occurrence of bloom depending on the nature of overall horizontal pull (triggered by k p , k h , u , v ) and vertical pull (triggered by w p , k v , w ) on phytoplankton biomass ( p 0 ) available on surface water, which are determined from the nature of dominating phytoplankton communities and variation in the nature of aquatic environment with seasonal changes.For turbulent medium (system 2), in the presence of high vertical mixing caused by higher vertical eddy diffusivity k v , when the vertical component of the spatial average of water velocity ( w ) dominates the corre- sponding horizontal component ( u , v ), being influenced by the change in the direction of the oceanic current, phytoplankton bloom of any community irrespective of whether it is buoyant ( k p > w p ) or its buoyancy gets dominated by its higher sinking tendency ( k p < w p ), can occur at a higher depth H c where sinking species gets stuck.Phytoplankton growth remains higher than phytoplankton loss at H c , that is, H c < D c , depending on the nature of vertical pull and grazing influence.If the dominating community tends to sink ( w p > k p , case BHD1), when w and vertical eddy diffusivity ( k v ) dominate over horizontal components ( u , v , k h ) in system 2, overall vertical pull on surface biomass increases.In such a scenario, dominating factors w p , k v , w take charge of the system's stability and the bloom's interconnected nature, provided that the grazing impact does not influence phytoplankton growth at the depth H c where sinking species accumulate.If the overall growth of phytoplankton communities at H c does not get annihilated by total loss over there, which happens when H c < D c , then the occurrence of bloom at H c depends only on the nature of dominating vertical velocity and diffusion parameters (Fig. 1(g)).
On the other hand, in the absence of a spatial average of water velocity, the chance of bloom formation at a higher depth ( H c ) for system 1 depends on the transportation rate of phytoplankton biomass along with the nutrient and the possibility of survival of sinking phytoplankton species at that depth ( H c ).In such a case, verti- cal eddy diffusivity ( k v ) plays a crucial role in bloom development at H c .
Stronger vertical pull influences the chances of bloom initiation at higher depths, but only higher values of vertical eddy diffusivity ( k v ), w , and w p are not sufficient enough for causing bloom at higher depth if mixing of biomass is limited to upper surface within the water column.This is possible in the presence of a high value of w , which induces a higher sinking tendency of dominating community, causing phytoplankton biomass to get distributed throughout the water column.However, the lack of vertical eddy diffusivity ( k v ) cannot provide the essential requirement of nutrients to the existing species at higher depths, which hinders species growth, implying that bloom cannot occur at higher depths.Higher vertical eddy diffusivity ( k v ) along with w being less than corresponding critical threshold values k v c , w c causes most of the phytoplankton biomass to get hetero- geneously and non uniformly distributed at different layers within the upper zone of the water column, which causes initiation of bloom at some of those layers, where grazing impact remains negligible and does not hinder phytoplankton growth.Figure 2c represents such a situation where low surface irradiance causes high vertical eddy diffusivity during the cold season in the presence of a high spatial average of water velocity along the vertical direction, caused by an exchange of oceanic current when higher sinking tendency of existing phytoplankton community dominates over its buoyancy ( w p > k p ).Current values of dominating factors ( k v , w ) being less than their critical values cause bloom ( H c , Fig. 2c ( 0 ≤ H c ≤ 5, 10 ≤ H c ≤ 20 )) or higher phytoplankton biomass (Fig. 2d ( 0 ≤ H c ≤ 5 )) at upper layers within the euphotic zone.One of these factors might be the reason for phytoplankton bloom development in the central Cantabrian Sea during late winter, where a slight increment in water temperature results in the reduction of vertical eddy diffusivity k v 14 .Besides, several phytoplankton patches of density varying from low (dark green) to high (yellow) observed in horizontal distribution (Fig. 2e) stand for the intermittency in the phytoplankton distribution of heavier communities during the cold season.This successfully explains observed intermittency in distributional dynamics of heavier species "Skeletonema Costatum" at a higher depth of 50 meters of Tokyo Bay in Dec 2006 and Feb 2008 8 .
Numerically, it has been observed that critical value k v c varies inversely with w irrespective of whether the community is buoyant ( k p > w p ) or has a higher sinking tendency ( w p > k p ) (Fig. 2(f)).Therefore, under the influence of the same circumstantial factors, the type of aquatic medium determined by the nature of water velocity decides the nature of phytoplankton bloom.Hence, in the presence of a high spatial average of water velocity along vertical direction satisfying w > w c for the turbulent environment (system 2), k v (> k v c ) can cause bloom at higher depth ( H c ) of water column (Fig. 1g, H c = L z = 60 meters), same k v in the absence of w , being extremely less than k v c in system 1, cannot transport excessive phytoplankton biomass from surface to deeper layers (Fig. 2d) and most of the phytoplankton biomass gets distributed within the upper surface of euphotic zones attached to the surface layer of water column (Fig. 2d), obstructing phytoplankton bloom initiation at higher depths.Therefore, if grazing by herbivorous zooplankton does not exert any negative impact on phytoplankton biomass available on surface water, then lack of transportation of excessive biomass from surface to bottom layer causes accumulation of phytoplankton biomass on the upper layer, which eventually generates the possibility of the occurrence of bloom at that upper layer during the winter-spring season for aquatic mediums like the lake, pond, and calm marine ecosystem, where wind velocity or exchange of seasonal currents are not the disturbing factors to ruin the calmness.On the other hand, if critical turbulence itself gets dominated by existing vertical eddy diffusivity ( k v > k v c ) and along with that, dominating component of the spatial average of water velocity remains higher than its critical value ( w > w c ) for turbulent medium (system 2), phytoplankton bloom can occur at a higher depth (say, H c ) or at the bottom layer of the water column ( L z ) (Fig. 1g), where intrinsic growth of phytoplankton is positive due to low grazing.
Though the bloom development at higher depth is triggered by higher vertical pull initially, it is not sufficient to cause bloom at deeper zones if the sinking rate of the existing phytoplankton community is negligible and the distribution of phytoplankton biomass is limited to the upper layers of the euphotic zone.This happens when existing species is almost buoyant ( w p ≈ 0 ), and it is the only dominating community during that season.Under such circumstances, higher vertical eddy diffusivity in the presence of a high spatial average of vertical water velocity causes transportation of some buoyant species into deeper levels for turbulent medium (system 2).However, in the absence of sufficient vertical pull, most phytoplankton biomass remains non-uniformly distributed at upper layers attached to surface water.Therefore, for system 2, even if the grazing impact remains low on surface water, such a non-intermittent assortment of existing biomass of dominating phytoplankton community within the upper layers of euphotic zone and sinking loss of existing community, both caused by vertical ) and the spatial average of water velocity ( w c ).When existing vertical eddy diffusivity k v > k v c along with w > w c for system 2, transportation of excessive mass to deeper level trig- gers chances of bloom at a higher depth H c or at the bottom layer L z (Fig. 1h), provided both H c or L z are less than critical depth D c .However, if k v < k v c and w < w c , then enough phytoplankton biomass of the buoyant community cannot be shifted to a deeper level to initiate bloom; instead, most of the phytoplankton biomass gets heterogeneously distributed within upper mixed layers and some sink to a specific small layer within the euphotic zone (Fig. 2g).Besides, this condition also triggers an upwelling event to transmit deep nutrient-rich water to the upper layer (Fig. 2h), which eventually causes nutrient abundance on the surface and upper oceanic layers generating the possibility of a phytoplankton bloom on surface water or in upper oceanic layers within the euphotic zone, where sinking phytoplankton biomass gets accumulated ( 0 ≤ H c ≤ 20 m as observed in Fig. 2g), and does not get negatively influenced by grazing pressure.This explains the fact of phytoplankton spring bloom initiation near-surface water of the temperate North Atlantic Ocean due to relaxation in turbulence profile 11 .It also supports the fact that phytoplankton bloom developed at a depth of nearly 10 meters of the mouth of the Arakawa river in May 2011 8 .
Critical vertical eddy diffusivity k v c increases as w decreases (Fig. 2f), as a result of which system 1 requires higher k v than that of the system 2 having non-zero w (see SI Tables 1, 2).Therefore, under the influence of the same circumstantial factors, existing vertical eddy diffusivity crossing its critical threshold value k v c transports sufficient biomass from the surface towards bottom layers in the presence of a high spatial average of vertical water velocity w greater than its critical value in system 2.However, the same vertical eddy diffusivity, in the absence of spatial average of vertical water velocity ( �w� = 0 ) for system 1, remains less than critical value k v c , which creates hindrance in bloom formation at higher depth due to the obstruction of such transportation of biomass (Fig. 3a).Besides, when it comes to the occurrence of bloom in a deeper zone of system 1, vertical eddy diffusivity (being higher than its critical value) triggers the chance of bloom initiation at a deeper level H c (Fig. 3b), provided H c is less than D c with low grazing.All the stated facts well explain the reason for phytoplank- ton bloom in higher depths of seas, oceans, or lakes during the cold season (Figs.1g,h, 3b,c), which has been observed in deep sea waters of North Arabian Sea during the winter season (Jan-March) 13 .
Numerical results under cases BHD1 and BHD2 also explain why surface bloom in winter is not so common in nature for system 2. First, lack of surface irradiance blocks phytoplankton growth on the surface mixed layer.
), (b) Vertical distribution of p 0 of buoyant community on t=40th day of cold season of system 1 ( k v > k v c , (c) Vertical distribution of p 0 on 95th day of cold season of system 1 ( k v > k v c ), p 0 has higher sinking tendency ( w p > k p ), (d) Vertical distribution of p 0 on 80th day of summer season of system 1 when p 0 is buoyant ( ), (e) Time variation of p 0 and z 0 on surface water ( k v << k v c ) of system 1 (summer season), p 0 is heavier ( w p > k p ), (f) Horizontal distribution of n 0 on surface layer on 62nd day of summer season ( k h > k h c , u > u c , v > v c ) for system 2, p 0 is buoyant in nature ( k p > w p ), (g) Horizontal distribution of n 0 on surface layer on 71th day of summer season ( k h > k h c , u < u c , v < v c ) of system 2, (h) Horizontal distribution of p 0 on 70th day of summer season, p 0 is buoyant in nature ( ) (red tide formation).The parameters values are mentioned in SI.
Besides, whatever the net growth is, extreme surface water stratification causes high vertical eddy diffusivity with k v > k v c when w > w c , influenced by the exchange of seasonal currents.This results in shifting most of the phytoplankton biomass to a deeper level (Fig. 1g,h), which ultimately blocks winter bloom on the surface mixed layer for most cases.As the winter season passes and early spring occurs, surface irradiance increases, which causes a slight increment in net productivity, resulting in a thicker phytoplankton layer on surface water that absorbs most of the incident solar irradiation.This explanation is valid for any aquatic environment regarding winter bloom formation at higher depths 8 .

Our next question is what triggers surface bloom formation for the almost calm Baltic Sea in spring season?
The study suggests that stratification of the water column is generally a prerequisite for most dinoflagellate blooms to develop in temperate areas 44,45 .In Baltic Sea, dinoflagellates are more sensitive to hydrographic conditions and climate fluctuation than to nutrient input in the Baltic sea on bloom onset 44 .Besides, experimental study has suggested that since the 1980s the nutrient loads to the Baltic Sea have been continuously decreasing 44 .According to this study 46 , an increase in thermal stratification can influence species-specific dinoflagellate distribution, behavior, and survival.Increased thermocline strength may decrease mixing between deep nutrient-rich and surface nutrient-depleted waters, leading to a decrease in surface water productivity 47 .However, stratification seems to favor species that can access the nutrient-rich water layers by the process of diel vertical migration (DVM), and the formation of algal blooms in low-nutrient surface waters is therefore possible 48 .Our study captures the dinoflagellates blooming on surface water in the spring season (Feb/March) in the southern Baltic Sea caused by water warming in the absence of extreme wind events 10 .The Baltic Sea is calm primarily during spring and summer, falling under system 1.As spring or summer arrives, water temperature increases due to an increment in irradiation, which reduces vertical eddy diffusivity k v .Such reduction causes k v to be less than its critical threshold value k v c .Hence, the transportation rate of mass from the surface to the bottom layer decreases for this almost calm marine environment in the absence of any wind event or exchange of oceanic current.Most phytoplankton species remain distributed on the upper layer within the water column.Therefore, species that have already sunk to higher depths and previously caused winter bloom at the depth H c (< D c ) die because of a lack of photosynthesis caused by a lack in the transportation of essential elements required for this process.Hence, the chance of blooming in higher depths reduces.In case the species is buoyant in nature, for example dinoflagellates, less sinking causes most of its biomass to accumulate on surface layer, which triggers the initiation of dinoflagellates surface bloom during early spring (Fig. 1a).In fact, excluding Baltic Sea, this explanation is valid for any aquatic medium under system 1 regarding the occurrence of surface bloom during early spring.Unless artificially induced vertical mixing generates instability in the water column and initiates the transportation of higher phytoplankton biomass and other ingredients from the surface to the bottom of the water column during the summer season, phytoplankton bloom usually occurs on the surface layer due to the impact of high thermal stratification of surface water.
Our third query, why is the surface bloom of the phytoplankton community so common in nature during summer, spring, or autumn?
The critical turbulence hypothesis proposes that a surface bloom can start under two different circumstances, (i) if the vertical eddy diffusivity is weak enough in the well-lit surface water for phytoplankton to receive sufficient light before being properly mixed beneath the critical depth or (ii) an intensely mixed surface layer shoals exposes phytoplankton to favorable light conditions 49,50 , which is possible when strong wind mixing period drives surface current to produce high horizontal eddy diffusivity ( k h ) in the presence of a high spatial average of the horizontal component of water velocity ( w ), which leads the existing community to get exposed to a favorable light and nutrient-rich zone.
As far as system 1 (lake, pond, or calm marine environment) is concerned, the absence of surface winds cannot trigger internal waves to inherit any horizontal turbulence.Hence, vertical eddy diffusivity is the only controlling factor for such a system, which acts on turbulent mixing and determines the nature of the planktonic distribution.When externally induced artificial mixing does not cause an imbalance in the nature of stationary flow, variation of vertical eddy diffusivity depends mostly on water density and temperature.During late spring, summer, and early autumn, high incident solar irradiance causes extreme thermal stratification to subjugate the density of water among water layers, resulting in reduced vertical eddy diffusivity k v .Under such circumstances, the nature of phytoplankton distribution varies depending only on the nature of critical vertical eddy diffusivity k v c .For system 1, zero contribution of the spatial average of water velocity to turbulent mixing requires k v c to be sufficiently higher for bloom initiation at higher depths, mostly when the buoyancy of existing community ( k p > w p ) drives most of the primary producers to float away and spread over a larger domain horizontally (Fig. 1(b)), than diving to deeper levels.Since k p > w p influences the horizontal pull, it stirs up the initiation of surface bloom in the presence of low grazing impact (Fig. 1(a)).When existing vertical eddy diffusivity k v (< k v c ) prevents the transportation of essential ingredients like temperature-induced warm water, nutrient, and several other abiotic factors required for photosynthesis, it accelerates the chance of initiation of surface bloom by causing no reduction in phytoplankton growth.This result explains why buoyant species float upwards and form dense surface blooms during weak mixing [51][52][53][54] .This result also mimics the fact that the buoyant phytoplankton community "Microcystis" developed a surface bloom during July 1992 in Lake Nieuwe Meer 15 .
Nevertheless, for stationary flow, phytoplankton bloom of buoyant class can also be initiated in an arbitrarily deeper layer (Fig. 3d) above critical depth in spring, early summer, or autumn season if externally induced artificial mixing triggers vertical eddy diffusivity to cross its critical threshold value.Then, existing surface biomass gets transported to a higher depth H c or to the bottom layer L z until its sinking to a higher depth of the water column is hindered by the water body's resilience.Then bloom can occur at that higher depth when k v > k v c , provided H c or L z < D c (Fig. 3d).This explains the fact of the occurrence of bloom of "Microcystis" at higher depth in Lake Nieuwe Meer during initial days in the presence of artificial mixing 15 .In fact, the hypothesis stated in 49 is supported by our numerical findings and the corresponding real-life occurrence regarding phytoplankton spring bloom initiation in an arbitrarily deep layer.
We now explain how k v < k v c controls the nature of the bloom of existing phytoplankton species if its sinking tendency dominates its buoyancy in system 1?
Fundamentally speaking, the heavier phytoplankton community sinks deeper, accelerating in the presence of higher vertical eddy diffusivity during the cold season.In general, strong vertical eddy diffusivity facilitates mixing, intensifying the possibility of bloom development of the heavier community at higher depth during that period (Fig. 1g).Though the higher sinking tendency of the heavier community like diatom triggers species sinking to deeper levels, the nature of vertical mixing sometimes controls and stops such sinking, depending on the nature of horizontal mixing and thermal stratification of water during the low wind mixing period in spring, early summer or autumn season.While wind speed does not take part in phytoplankton distribution during the thermally stratified period at the beginning of the spring season in lake water (system 1), increment in water temperature hinders increment in the dominating factor k v while the spatial average of the vertical velocity of water is zero.This results in k v < k v c , which obstructs the transportation of surface biomass, including nutrients, several organic factors, and phytoplankton biomass from the surface to the bottom layer.Besides, the occurrence of vertical upwelling event influenced by the fact k v < k v c , causes vertical transmission of nutrients and sinking phytoplankton biomass to the surface layer, as a result of which deep nutrient-rich water reaches the surface level, which provides an ample supply of nutrient to the existing phytoplankton community that enhances phytoplankton productivity and growth.An abundance of nutrients (Fig. 3e) and high water temperature accelerate phytoplankton productivity and growth.Low wind mixing generates low horizontal eddy diffusivity k h , which causes biomass to get distributed over surrounding regions and reduces the pressure of phytoplankton biomass on the surface layer.Hence, high productivity and spreading of phytoplankton biomass initiates surface bloom development (Fig. 1c).Under such circumstances, phytoplankton bloom gets initiated for a specific time, even when the impact of high grazing pressure causes zooplankton abundance (Fig. 2b).In fact, this result mimics the event of the occurrence of bloom of both phytoplankton class (diatom with higher sinking tendency in shallower lakes) and zooplankton class (Cladocerans dominated by the Daphnia hyalina-galeata complex) as mentioned in 16 .

Our next analysis answers the question what drives surface bloom formation of buoyant phytoplankton community in a turbulent marine environment?
Wind-induced horizontal turbulence is one of those critical factors for developing surface bloom in system 2. In a turbulent environment, wind-driven horizontal currents deliver deep nutrient-rich water to the surface photic layers.Also, in the presence of weak vertical eddy diffusivity during the warmer season, wind-induced strong horizontal eddy diffusivity ( k h ) triggers horizontal mixing, which along with a higher spatial average of horizontal water velocity ( u , v ) of system 2 lead all plankton bodies, nutrients, warmer water to a compara- tively stationary surface zone.A fair wind is required to cause such transportation since only higher values of horizontal water velocity and turbulence are not sufficient enough to initiate surface bloom.Corresponding to that, we have considered there exist critical values k h c , u c , v c (say) of horizontal eddy diffusivity and the spatial average of horizontal water velocity component, respectively, which decide whether the existing surface wind is appropriate to trigger the surface bloom.If existing horizontal eddy diffusivity dominates its critical value, while the spatial average of horizontal water velocities have already crossed their critical values, horizontal transportation of phytoplankton biomass along with required ingredients from dense, turbulent zone to stationary zone ) generates the possibility of surface bloom initiation (Fig. 1e) in the presence of low grazing pressure on surface water.
On the contrary, whenever this condition is violated by the nature of wind ( ), lack of horizontal transportation of surface biomass makes it heavier to dominate water resilience and sink into slightly deeper zone getting influenced by existing vertical turbulence.However, weak vertical mixing caused by weak vertical eddy diffusivity ( k v < k v c ) during such a thermally stratified period with w < w c cannot drive phytoplankton species deeper; instead, they sink to a particular layer ( H c ) within the upper zone of the euphotic layer where nutrients and light are sufficiently available for the continuation of photosynthesis to maintain marine food web.This results in higher phytoplankton biomass at that layer (Fig. 1d, H c = 20 m) with negligible grazing.
Generally, the effect of turbulent current is different for different phytoplankton species.In particular, the conditions stand for the bloom of those species which are buoyant and capable of growing in a highly turbulent environment.Wind-induced horizontal turbulence drives nutrients along with all patches of buoyant species to an exceptionally stable surface layer, where the currents are more or less steady, and vertical turbulence is critically weak due to extreme thermal stratification triggered by intense solar irradiation during summer, late spring or early autumn season and nutrient density keeps on growing on the surface layer with time (Fig .3f→ g).Therefore, when grazing pressure does not exert any negative impact on phytoplankton on surface water (Fig. 2a), the surface bloom of the existing buoyant community gets developed in nutrient-rich saline water (Fig. 1e).This might possibly be the reason behind several nonmotile flagellate classes, like "Alexandrium Tamarense", forming red tides (Fig. 3h) (flagellates bloom) in Mikawa Bay of Japan in a severe diffusive condition with horizontal turbulence k h is of order 10 2 m 2 s −1 and vertical turbulence k v is of order 10 −5 m 2 s −1 in spring, 1991 (28th March-22nd April) 9 .Whenever the growth of turbulent driven species is immensely responsive to such an environment ( ), surface bloom or red tides (Fig. 3h) of motile buoyant classes gets formed as well in the presence of low grazing impact.Such a scenario explains the red tide formation of several motile, buoyant flagellates like "Ceratium Furca" (July 1990: Summer), "Ciliate", "Mesodinium Rubrum" (November 1991: Autumn), etc. in Mikawa Bay, Japan during several periods from April 1989 to Dec 1991 9 .Now, we would like to address what happens to the occurrence of surface bloom while dominating class has a higher sinking tendency (w p > k p ) due to heavier size in a wind-induced turbulent flow?
w p triggers vertical pull on available biomass of heavier community; hence, ensuring the possibility of sur- face bloom development requires intense wind-driven horizontal eddy diffusivity, which will transport heavier www.nature.com/scientificreports/phytoplankton body to a dense nutrient-rich zone to dominate its sinking by the resilience of the dense water body of that zone.Generally, during a thermally stratified period, vertical eddy diffusivity reduces and reaches beneath its critical level ( k v < k v c ), as a result of which existing k v cannot transport excessive mass from surface to bottom layer; instead, it initiates upwelling of water, which transports nutrient-rich water of bottom layers to surface level.As this happens, an increment in water density enhances the chance of water column instability, which triggers two related events under two different circumstances.First, during intense wind mixing period, if existing horizontal eddy diffusivity dominates its critical value ( k h > k h c ) along with u > u c , v > v c ), then this drives extreme dense water along with existing phytoplankton biomass from highly dense, turbulent zone to approximately stable region, where dense water and overall biomass spread over a larger domain to maintain system stability.As the overall biomass of nutrients and phytoplankton gets assembled in a stable domain, the resilience of dense water bodies eventually obstructs species from sinking to higher depths.Besides, in the presence of low grazing impact, increment in nutrient biomass with time provides an ample supply of nutrients (Fig. 4a→ b), here the light yellowish shaded zone stands for higher nutrient biomass in Fig. 4b compared to Fig. 4a, to the existing phytoplankton community on the surface zone.Additionally, inadequate light availability on the surface layer triggers species productivity which eventually strengthens the possibility of surface bloom (Fig. 5c,f).Such bloom of heavier species, namely, dinoflagellates and diatoms, has been observed in the southern Baltic Sea in Feb/March 10 .The above facts also explain the event of spring bloom initiation on the surface water of the North Atlantic Ocean during a thermally stratified period associated with strong westerly winds in winter-spring 55 .
It is highly unrealistic for a wind-driven turbulent medium (system 2) in the absence of any intense exchange of oceanic currents (implying �w� ≈ 0 ) to have bloom at higher depth, since during wind mixed and stratified period k h remains of order 10 2 , therefore, the condition k v > k v c to hold, k v needs to be of order 10 15 to cross k v c as it has numerically been observed for k h = 9 × 10 3 , k v c ≥ 2.3658 × 10 14 , which is practically impossible as thermal stratification does not let k v to increase this much unless the system is almost stationary where k h will be critically low.Therefore, for most cases, even if k h < k h c causes an increment in k v , it remains less than k v c unless ) of system 2, p 0 is heavier in nature ( k p < w p ), (b) Horizontal distribution of n 0 on surface layer on 70th day of summer season with same conditions, (c) Horizontal distribution of phytoplankton variance �p ′2 � on 75th day of warmer season of system 2 ( ), p 0 is heavier ( w p > k p ), (d) Time variation of p 0 and z 0 in summer season of system 2 ( k v < k v c ), p 0 is heavier, (e) Vertical distribution of p 0 of heavier community on 168th day of summer season of system 2 ( ), (f) Vertical distribution of z 0 on 168th day of summer season ( I 0 = 220 W/m 2 ) of system 2, with same condition.All the figures are in the presence of high grazing influence.The parameters values are mentioned in SI.
the turbulency is gone due to eventual stopping in wind event during a thermally stratified period in the absence of any intense exchange of oceanic currents.Whenever this happens, though heavier weight causes species to dive deeper, the fact k v < k v c along with w < w c hinders transportation of surface biomass to a deeper level, and whenever some species sink due to sinking nature, they get shifted to the upper layer and remains heterogeneously distributed over there (Fig. 1f).In particular, variation in the nature of the horizontal eddy diffusivity and water velocity generates the possibility of two different circumstances, (i) , which can be biologically interpreted as the impact of the episodic wind event.When wind hits the ocean surface, sudden acceleration in the horizontal eddy diffusivity and water velocity results in , which initiates a bloom cycle.However, in between two wind events, a particular gap exists that causes a reduction in horizontal eddy diffusivity and velocity distribution of oceanic flows.Such reduction causes k h < k h c , u < u c , v < v c , which makes the system unstable during the bloomy period, and species starts to sink again.Under such circumstances, depending on the nature of k v , w and driven by the fact whether ( k v < or > k v c , w < or > w c ), species remains homogeneously distributed in the upper ocean due to upwelling of water ( k v < k v c , w < w c (Fig. 1d,f)) or they dive deeper due to k v > k v c , w > w c (Fig. 5a,b).If another episodic wind event hits the ocean surface during this process, then it again activates horizontal components to cause k h > k h c , u > u c , v > v c and k v < k v c , w < w c .This breaks the cycle of species sinking, and the bloom cycle begins again similarly (Fig. 5d-f).
This bloom cycle has been numerically captured in Fig. 5, where phytoplankton bloom has been observed on surface water and the upper layer of the euphotic zone on the 62nd and 75th day of the summer season in two weeks.For the first bloom, starting from the 50th day onwards, the sinking of heavier species to deeper levels eventually starts to decrease, and for the second bloom, starting from the 65th day onwards, the sinking of heavier community decreases, which happens due to variation in the nutrient density on surface water in the presence of episodic wind event.This is observed numerically in the case SB4 where nutrient density on the surface layer is higher on the 70th day (Fig. 4b) than that on the 64th day (Fig. 4a), which makes water nutrient-rich in the presence of low grazing. (a) Vertical distribution of p 0 on 50th day of cold season of system 2 ( k h < k h c , u < u c , v < v c ), (b) Vertical distribution of p 0 on 53rd day of cold season with same conditions, (c) Vertical distribution of p 0 on 62nd day of cold season of system 2 ( k h > k h c , u > u c , v > v c ), (d) Vertical distribution of p 0 on 65th day of cold season of system 2 ( k h < k h c , u < u c , v < v c ), (e) Vertical distribution of p 0 on 70th day of cold season of system 2 with same conditions, (f) Vertical distribution of p 0 on 75th day of cold season of system 2 ( ).The parameters values are mentioned in SI. and provides ample nutrients required for photosynthesis to the existing phytoplankton community on surface layer resulting enhancement in phytoplankton community growth.Generally, the sinking rate of heavier classes like diatom increases with nutrient stress 56 .As a result, an inadequate supply of nutrients driven by the fact k h > k h c , u > u c , v > v c reduces nutrient stress and hinders sinking of such heavier class, which entertain the possibility of bloom formation of heavier community on surface layer (Fig. 5c).This bloom now sustains for three days, and again due to the eventual stopping of wind event ( k h < k h c , u < u c , v < v c ), excessive phytoplankton biomass cannot spread over the larger domain and triggers water column instability.As a result, in the absence of a wind event, when there is no intense exchange of oceanic current ( �w� ≈ 0 ), k v increases to restore water column stability but thermal stratification cannot let k v cross k v c , therefore on the very first day, species get transported to deeper layer (65th day, (Fig. 5d)), but eventually k v < k v c , w < w c triggers upwelling event, as a result of which amount of sinking biomass reduces (Fig. 5e), as wind strikes again on 70th day, the fact k h > k h c , u > u c , v > v c keeps on providing ample supply of nutrient to existing phytoplankton com- munity, and phytoplankton productivity gets enhanced.Under such conditions, extreme phytoplankton biomass floats away towards a more extensive domain (Fig. 4c), and water column stability, even in the presence of higher phytoplankton biomass, remains intact.Eventually, bloom occurs on the upper ocean surface (Fig. 5f).This also explains the fact why the bloom cycle of heavier species "Skeletonema Costatum" has been observed within a gap of certain days during wind mixing and thermally stratified period of 1990, 1991 (11-20  We conclude the discussion by answering our final question, how does grazing control surface bloom dynamics of phytoplankton on two opposite sides of critical vertical eddy diffusivity k v c for the almost calm aquatic environment?Grazing by herbivorous zooplankton plays a vital role in both scenarios, whether the bloom occurs on the surface layer or at a higher depth.For a fixed region with seasonal changes, the breeding and reproduction of herbivorous zooplankton generate differences in the overall population density of the grazing community, which controls the overall density of the phytoplankton community.In fact, at a particular time, this situation might arise with regional differences.As discussed in case SB3, whenever critical horizontal eddy diffusivity gets dominated by existing wind-induced horizontal eddy diffusivity in a turbulent flow, where the spatial average of intense horizontal velocity dominates its critical value, the chance of surface bloom development increases.Surface bloom under such circumstances can only be developed if and only if grazing impact remains negligible (as observed in SB3, SB4, Fig. 1e, 5c,f).In case, the grazing rate and assimilated grazing coefficient in the presence of low intra-specific competition among zooplankton classes exert a negative impact on the phytoplankton growth curve (Fig. 4d), formation of phytoplankton bloom (Fig. 4e) on surface water gets obstructed by higher zooplankton grazing on surface layer (Fig. 4f), even in a compatible environment in terms of horizontal eddy diffusivity and the spatial average of water velocity w .It should be noted that higher grazing pressure cannot successfully dominate phytoplankton bloom development on the surface layer for a certain period for those aquatic environments, where an ample amount of nutrient gets supplied to the existing phytoplankton community in the presence of a ongoing upwelling event triggered by the condition k v < k v c , w < w c for system 2 and k v < k v c for system 1, which not only supplies sufficient nutrient but also transmit sinking phytoplankton biomass to surface layer.An abundance of nutrients and lack of sinking of phytoplankton class generates surface bloom even in the presence of higher grazing pressure for a certain amount of time, which has been observed in case SB2 (Fig. 1c), when compared to zooplankton biomass, nutrient density is so high that phytoplankton and zooplankton both grow together between approximately 70-110 days of the warmer season (Fig. 3e).Therefore, phytoplankton bloom gets developed on the surface layer (Fig. 1c) on the 115th day, while zooplankton is also abundant on the same day on the surface layer (Fig. 2b) due to high grazing pressure.Hence, if the nutrient supply remains much higher, then there exists a possibility of phytoplankton bloom for a certain amount of time, but, in case nutrient supply and zooplankton biomass both remain high together in a nutrient-rich zone during no wind mixing period and in the absence of any intense exchange of oceanic current ( �w� ≈ 0 ), an absence of upwelling event does not transport deep nutrient-rich water and sinking phytoplankton biomass to the upper surface, instead k v > k v c shifts excessive phytoplankton biomass of buoyant community to a deeper layer.Therefore, in one direction, ample nutrient supply in nutrient-rich zones causes richness in phytoplankton biomass even in the presence of higher grazing pressure.On the other hand, as soon as excessive growth disturbs water column instability, surface biomass gets transported to the deeper layer.Whatever phytoplankton remains on the surface layer, higher grazing hinders its bloom formation on surface water (as observed in Fig. 4e).On the contrary, the availability of inadequate food sources in higher grazing pressure results in zooplankton bloom development on surface water (Fig. 4f).This has been numerically captured in Fig. 4, where k v > k v c along with high grazing pressure due to the abundance of zooplankton on the surface layer (Fig. 4f) obstruct phytoplankton growth (Fig. 4e) on surface water in a zone where both zooplankton and nutrient biomass are higher together (Fig. 4d).In fact, such a situation has been observed in a high nutrient zone of Gulf Alaska, where low diatom density is followed by an abundance in zooplankton 12 .

Conclusion
The occurrence of phytoplankton bloom has always been a global phenomenon throughout all seasons for most aquatic ecosystems, including marine life of sea or ocean or stationary systems like lakes or ponds.Spatial and temporal contrasts in temperature, salinity, and chemical composition, along with differences in the distribution of biological agents, characterize the nature of aquatic systems.Depending on the quality of the water body and the characteristic of aquatic life, the nature of the phytoplankton community varies regionally and seasonally.Generally, there is a critical threshold value of water density for any aquatic system, whether it is freshwater or marine water.When this value is crossed, water becomes extremely stratified causing an imbalance in vertical mixing.This triggers vertical turbulent diffusivity k v and in the long run affects the distribution of the phytoplank- ton community.Habitually, variation in water temperature due to the lack of thermal stratification caused by low surface irradiance, wind-induced evaporation, and snowfall causes variation in water density.For example, the nature of thermal stratification governed by surface irradiance fluctuates with atmospheric variation and water density on the surface layer, which varies inversely with this thermal stratification and is directly proportional to salinity.The quality of water density stipulates the level of stratification.Generally, water quality varies with seasonal changes; hence level of stratification varies with climate change.As this happens, the stability of the water column gets affected, which influences the nature of horizontal ( k h ) and vertical eddy diffusivity ( k v ).This vari- ation in diffusivity causes variation in mixing in a turbulent flow.On the other hand, changes in the direction of oceanic current with seasonal changes causes disturbances in water temperature.The surface water temperature remains relatively low during winter due to less thermal stratification, which causes an increment in the density of surface water.When this density crosses a critical threshold value, higher density on the surface layer makes the water column unstable.As a result, vertical turbulent mixing increases to transport excessive mass of plankton along with nutrients, which is required for photosynthesis, carried within the water column from the upper layer to the bottom layer until the column gets stabilized (Figs.1g,h, 2d,h).It has long been recognized that the conditions of turbulence and mixing are critical factors for the growth and persistence of natural populations of phytoplankton in lakes and oceans [57][58][59] .
It should be noted that excluding eddy diffusivity k v , k h , the nature of water advection, being influenced by the spatial average of water velocity, also plays a vital role in the distribution of phytoplankton biomass.As far as systems like the lake, pond and calm oceanic medium (system 1) are concerned, usually spatial average of horizontal ( u , v ) and vertical ( w ) water velocities remain zero in the absence of any externally induced artificial mixing.However, for a turbulent medium like a flowing river or rough marine ecosystem (system 2), external forces like seasonal tides, wind velocity, oceanic current, and sometimes natural calamity cause disturbances in the horizontal and vertical velocity distribution of water, hence the spatial average of horizontal and vertical water velocity, �u� � = 0, �v� � = 0, �w� � = 0 , as a result of which distribution of planktonic system gets influenced.
The distribution of phytoplankton biomass and chances of bloom occurrence depend mainly on dominating nature of the control parameters k v , k h , k p , w p , u , v , and w when the grazing impact is low.Among these control parameters, horizontal components ( k h , k p , u , v ) trigger the chances of surface bloom, whereas vertical components ( k v , w p , w ) enhance the chance of having bloom at higher depth.Dominating factors have specific critical values (determined from Routh Hurwitz stability criteria).The nature of the bloom and stability of the system depends on whether the dominating factors have crossed their corresponding critical threshold value in particular circumstances determined by the parameter values of the system.If each of the dominating factors is higher than the corresponding critical value, this determines the nature of phytoplankton bloom.However, once the bloom is developed either on the surface layer or at a higher depth, its sustainability depends on the system's stability in particular circumstances determined by the parameter values of the system.As soon as the stability of the system is assured, if phytoplankton biomass remains higher at a certain depth H c or at bottom layer L z satisfying H c or L z < D c under the considered circumstances, then this higher biomass will trigger bloom initiation at that depth H c or at the bottom layer L z of shallow water being influenced by the stability of the system, which has been observed for all above considered cases (SB1-SB4, BHD1, and BHD2).
Therefore, for any aquatic environment like a still lake, pond, or almost calm marine environment (system 1), we may conclude the following facts regarding phytoplankton bloom dynamics in the presence of low grazing impact, (i) surface bloom of buoyant species will mostly occur during a warmer season without any artificial or wind-induced mixing.Even in the presence of artificial mixing, a surface bloom of buoyant species can occur when vertical eddy diffusion k v < k v c cannot transport all phytoplankton biomass to higher depth (Case SB1, (SI Table 1)).This situation is observed when multi-species exist in a lake, pond, or almost calm marine environment, where one community is buoyant, and some have a higher sinking tendency.Then, during the warmer season, there will be a higher chance of having surface bloom of buoyant species under the circumstances where existing vertical eddy diffusion remains sufficiently low to be dominated by its corresponding critical turbulence.(ii) But, if artificial mixing or cold weather influences the aquatic environment of the lake, pond, or almost stationary aquatic flow and affects the stability of the water column, then most bloom will occur at higher depths in the presence of low grazing depending on two facts (a) dominance of k v c by existing vertical eddy diffusivity k v (b) the depth ( H c or L z ) at which sinking species gets stuck, is less than the critical depth ( D c ) (cases BHD1, BHD2, SB1, (SI Table 1).
As far as turbulent flow is concerned (system 2), we obtain the following facts regarding bloom formation in the presence of low grazing influence, (i) when wind-induced horizontal turbulent diffusion k h dominates vertical turbulence k v mostly during a thermally stratified period and the spatial average of horizontal water velocity dominates corresponding vertical component as well influenced by the act of surface wind event.The chance of surface development, irrespective of the nature of the existing community, increases provided the fact that horizontal eddy diffusivity and the spatial average of horizontal water velocity dominate their corresponding critical values ( k h > k h c , u > u c , v > v c , cases SB3, SB4, 1e, 5d,f (SI Table 2)).If the system is in motion, and the horizontal mixing rate dominates, then the surface bloom of any species can also occur in turbulent k p )P − ∂ ∂y (v + k p )P − ∂ ∂z (w + w p )P + S 2 (N, P, Z),

Figure 1 .
Figure 1.The figure shows the distribution of phytoplankton biomass ( p 0 ) under different circumstances.(a)Vertical distribution of p 0 of buoyant community of system 1 ( k v < k v c ), (b) Horizontal distribution of p 0 of the buoyant community of system 1 (k v < k v c ), (c) Vertical distribution of p 0 of system 1 ( k v < k v c ), (d) Vertical distribution of p 0 of buoyant community of system 2 ( k v < k h , k h < k h c , u < u c , v < v c ), (e) Vertical distribution of p 0 of buoyant community of system 2 ( k v < k h , k h > k h c , u > u c , v > v c , (f) Vertical distribution of p 0 of heavier community of system 2 ( k v < k h , k h > k h c , u > u c , v > v c, (g) Vertical distribution of p 0 of heavier community of system 2 ( k v > k v c , w > w c ) (h) Vertical distribution of p 0 of buoyant community of system 2 ( k v > k v c , w > w c ).All the figures are in the presence of low grazing influence.The parameters values are mentioned in SI.

Figure 2 .
Figure 2. (a) Time variation of p 0 and z 0 in summer season of system 2, p 0 is buoyant in nature( k p > w p ), k h > k h c , u > u c , v > v c , (b)Vertical distribution of z 0 on 115th day of summer season of system 1, p 0 is heavier (w p > k p ), ( k v < k v c), (c) Vertical distribution of p 0 on 60th day of cold season ( k v < k v c , < w >< w c ) of system 2, p 0 has higher sinking tendency ( w p > k p ), (d) Vertical distribution of p 0 on 48th day of cold season of system 1 ( k v << k v c ), p 0 has higher sinking tendency, (e) Horizontal distribution of p 0 at the depth H c = 60 meters on 65th day of cold season ( k v > k v c , < w >> w c ) of system 2, p 0 has higher sinking tendency, (f) Relation between w and k v c using estimated data set, (g) Vertical distribution of p 0 of buoyant community ( k p > w p ) on 35th day of cold season of system 2 ( k v < k v c , w < w c ), (h) Vertical distribution of n 0 of system 2 on 32nd day of cold season ( k v < k v c , w < w c ), p 0 is buoyant in nature ( k p > w p ).The parameters values are mentioned in SI. https://doi.org/10.1038/s41598-023-43745-z https://doi.org/10.1038/s41598-023-43745-zwww.nature.com/scientificreports/

Figure 3 .
Figure 3. (a) Vertical distribution of p 0 of buoyant community on t=35th day of cold season of system 1( k v < k v c), (b) Vertical distribution of p 0 of buoyant community on t=40th day of cold season of system 1 ( k v > k v c , (c) Vertical distribution of p 0 on 95th day of cold season of system 1 ( k v > k v c ), p 0 has higher sinking tendency ( w p > k p ), (d) Vertical distribution of p 0 on 80th day of summer season of system 1 when p 0 is buoyant (k p > w p , k v > k v c), (e) Time variation of p 0 and z 0 on surface water ( k v << k v c ) of system 1 (summer season), p 0 is heavier ( w p > k p ), (f) Horizontal distribution of n 0 on surface layer on 62nd day of summer season ( k h > k h c , u > u c , v > v c ) for system 2, p 0 is buoyant in nature ( k p > w p ), (g) Horizontal distribution of n 0 on surface layer on 71th day of summer season ( k h > k h c , u < u c , v < v c ) of system 2, (h) Horizontal distribution of p 0 on 70th day of summer season, p 0 is buoyant in nature ( k p > w p , k h > k h c , u > u c , v > v c ) (red tide formation).The parameters values are mentioned in SI. https://doi.org/10.1038/s41598-023-43745-z

Figure 4 .
Figure 4. (a) Horizontal distribution of n 0 on surface layer on 64th day of summer season( k h > k h c , u > u c , v > v c) of system 2, p 0 is heavier in nature ( k p < w p ), (b) Horizontal distribution of n 0 on surface layer on 70th day of summer season with same conditions, (c) Horizontal distribution of phytoplankton variance �p ′2 � on 75th day of warmer season of system 2 (k h > k h c , u > u c , v > v c), p 0 is heavier ( w p > k p ), (d) Time variation of p 0 and z 0 in summer season of system 2 ( k v < k v c ), p 0 is heavier, (e) Vertical distribution of p 0 of heavier community on 168th day of summer season of system 2 (k h > k h c , u > u c , v v c , k v > k v c), (f) Vertical distribution of z 0 on 168th day of summer season ( I 0 = 220 W/m 2 ) of system 2, with same condition.All the figures are in the presence of high grazing influence.The parameters values are mentioned in SI.

Figure 5 .
Figure 5. Bloom cycle of heavier community of p 0 as k h< k h c , u < u c , v < v c shifts to k h > k h c , u > u c , v > v cin the presence of low grazing. (a) Vertical distribution of p 0 on 50th day of cold season of system 2 (k h < k h c , u < u c , v < v c), (b) Vertical distribution of p 0 on 53rd day of cold season with same conditions, (c) Vertical distribution of p 0 on 62nd day of cold season of system 2 ( k h > k h c , u > u c , v > v c ), (d) Vertical distribution of p 0 on 65th day of cold season of system 2 (k h < k h c , u < u c , v < v c), (e) Vertical distribution of p 0 on 70th day of cold season of system 2 with same conditions, (f) Vertical distribution of p 0 on 75th day of cold season of system 2 ( k h > k h c , u > u c , v > v c ).The parameters values are mentioned in SI.

Table 2 .
Dimension and ranges of different quantities used in this model and their dimensions.