A two-timescale model of plankton–oxygen dynamics predicts formation of oxygen minimum zones and global anoxia

Decline of the dissolved oxygen in the ocean is a growing concern, as it may eventually lead to global anoxia, an elevated mortality of marine fauna and even a mass extinction. Deoxygenation of the ocean often results in the formation of oxygen minimum zones (OMZ): large domains where the abundance of oxygen is much lower than that in the surrounding ocean environment. Factors and processes resulting in the OMZ formation remain controversial. We consider a conceptual model of coupled plankton–oxygen dynamics that, apart from the plankton growth and the oxygen production by phytoplankton, also accounts for the difference in the timescales for phyto- and zooplankton (making it a “slow-fast system”) and for the implicit effect of upper trophic levels resulting in density dependent (nonlinear) zooplankton mortality. The model is investigated using a combination of analytical techniques and numerical simulations. The slow-fast system is decomposed into its slow and fast subsystems. The critical manifold of the slow-fast system and its stability is then studied by analyzing the bifurcation structure of the fast subsystem. We obtain the canard cycles of the slow-fast system for a range of parameter values. However, the system does not allow for persistent relaxation oscillations; instead, the blowup of the canard cycle results in plankton extinction and oxygen depletion. For the spatially explicit model, the earlier works in this direction did not take into account the density dependent mortality rate of the zooplankton, and thus could exhibit Turing pattern. However, the inclusion of the density dependent mortality into the system can lead to stationary Turing patterns. The dynamics of the system is then studied near the Turing bifurcation threshold. We further consider the effect of the self-movement of the zooplankton along with the turbulent mixing. We show that an initial non-uniform perturbation can lead to the formation of an OMZ, which then grows in size and spreads over space. For a sufficiently large timescale separation, the spread of the OMZ can result in global anoxia.


Introduction
Plankton is a vital element in the complex marine food webs and biochemical cycles.Phytoplankton is the primary producer standing at the base of the marine food web.As a by-product of the primary production, phytoplankton produces oxygen in the process called photosynthesis.The produced oxygen is then used by marine fauna, e.g.zooplankton and fish.The level of dissolved oxygen is a crucial indicator of the marine ecosystem health, as its depletion may lead to a mass mortality of aquatic species [1][2][3].Furthermore, a considerable part of the oxygen produced in the ocean goes to the atmosphere through the ocean surface.It is estimated that around 50-80% of atmospheric oxygen originates in the ocean.Phytoplankton therefore plays a critical role in producing and maintaining oxygen levels needed for survival both of aquatic and terrestrial species [4,5].
Over the last few decades, the level of the dissolved oxygen in the ocean has shown a trend to decrease [6][7][8][9][10].This is believed to be a result of the global warming, in particular because warmer water contains less dissolved oxygen [6,11].Also, the warming leads to a stronger stratification of the upper ocean, which reduces the O 2 fluxes through the ocean surface [12,13].Apart from the above purely physical mechanisms, there can be more subtle effects of the global warming driven by a biological feedback.An increase in the water temperature may slow down the oxygen production by phytoplankton [14][15][16][17], potentially resulting in a regime shift and a global oxygen depletion [18][19][20].Ocean anoxia is thought to be the factor that can trigger a mass extinction and this has indeed happened several times during the deep past [3,[21][22][23].Thus, better understanding of the pathways leading to the global anoxia as well as the identification of possible early signs of the approaching catastrophe are obviously problems of literally vital importance.In turn, it requires a better understanding of the coupled plankton-oxygen dynamics in the ocean and, arguably, mathematical modelling is a powerful research approach to facilitate it.
Mathematical models of plankton dynamics are abundant in the literature, e.g.see [24][25][26][27][28][29][30][31].Earlier modelling studies of the plankton dynamics were mostly concerned with the temporal and spatio-temporal dynamics of the coupled phytoplankton-zooplankton system [28,[32][33][34], with a particular focus on plankton patchiness and plankton blooming [24,31,[35][36][37].In [38][39][40][41][42], conceptual two-and three-component mathematical models were considered to reveal the role of various internal and external factors and to demonstrate different routes to plankton pattern formation.More recently, there has been growing attention to possible links between the dissolved oxygen, plankton dynamics and the climate change [43,44], which facilitated further progress in mathematical modelling of marine ecosystems [45].In particular, it was shown in [18,19,46] that sustainable oxygen production by marine phytoplankton can be severely disrupted by a gradual increase in the average water temperature.In turn, this may eventually lead to a global anoxia and it was argued in [47] that the observed decrease in the oxygen stock in the ocean [6,7,48] and the slow gradual decrease in the amount of atmospheric oxygen [49] that has occurred over the last few decades may be early signs of the approaching catastrophe.
Remarkably, the amount of dissolved oxygen does not only change with time, it also depends on space.The spatial distribution of oxygen in the ocean is distinctly heterogeneous [50,51], sometimes resulting in the formation of large stable areas or 'patches' where the dissolved oxygen concentration is much lower than the average.Such a patch is referred to as Oxygen Minimum Zone (OMZ) or the dead zone [2,7,52].They were discovered in different parts of the world ocean, e.g. in the subsurface waters of the Arabian Sea and in the eastern boundary upwelling regions of the tropical oceans of California, Peru, and Namibia [52,53].The existence of OMZ has a significant effect on the marine species abundance and the aquatic food chains [54].There is evidence that some zooplankton species may have a capacity to adapt to oxygen-deficient environments [55].However, a significant drop in the dissolved oxygen level eventually results in the formation of a dead zone, so that the majority of marine life either dies or leaves the area [2,54].
Interestingly, over the last several decades the OMZs have been growing in size.In particular, a rapid expansion in OMZs in the eastern Pacific and northern Indian oceans is well documented [7,52].Because of the OMZ's detrimental effect on the corresponding local marine ecosystem, this is becoming a grave concern for ecology and conservation as well as some industries, e.g.fishery.Moreover, there is some theoretical evidence that the OMZ expansion may be an early warning signal of the approaching global anoxia [56].The expansion of OMZ is thought to be caused by various factors, ultimately linking it to the global warming and to the human interference through the perturbation of the ocean's biogeochemical cycles, in particular the CO 2 cycle, although the issue as a whole remains controversial [8,57].
In this study, we address the phenomena of the OMZ formation and growth theoretically based on the earlier conceptual modelling approach that considers the variations (in particular, a decrease) in the dissolved oxygen level as a inherent property of the coupled plankton-oxygen dynamics in the ocean, not necessarily an effect of exogenous factors [18][19][20].Our updated mathematical model incorporates two important factors that were overlooked in the earlier studies.One such factor is the nonlinear mortality of zooplankton; it takes into account a combined effect of the zooplankton intraspecific competition, cannibalism, and the zooplankton consumption by its predators from a higher trophic level (e.g.fish) [28,38].It is well known that the nonlinear mortality rate can change the system's dynamical properties significantly [31,58,59] but its potential effect on the plankton-oxygen dynamics has never been investigated.
The second factor is the existence of different timescales in the plankton-oxygen dynamics.Indeed, it is a common observation that the zooplankton growth rate is usually much lower than that of phytoplankton.Correspondingly, the typical time (timescale) of changes in the phytoplankton density is considerably shorter (sometimes by an order of magnitude) than that of zooplankton.In this study, we therefore assume that the production of oxygen and the phytoplankton growth, on the one hand, and the zooplankton response, on the other hand, happen on fast and slow timescales, respectively.The presence of different timescales in a dynamical system can make its properties much more complicated, e.g. to bring bifurcations, coexistence of multiple attractors, complex oscillations, long transients and pattern formation that would not be there otherwise [60][61][62][63][64].It is therefore a very relevant question as to how the existence of multiple timescales can modify the plankton-oxygen dynamics, in particular in the context of the OMZ formation and growth.
The paper is organised as follows.In the next section, we describe our mathematical model and investigate its basic properties such as the existence and stability of the steady states.In Sections 3 and 4, we consider the properties of the nonspatial model to reveal, respectively, the effect of the nonlinear mortality and the different timescales.We identify conditions when the system can undergo a regime shift that may correspond to catastrophic changes in the real world.In Section 5, we consider the properties of the spatially explicit system, with a special focus on the dynamical regimes resulting in the pattern formation, in particular those that can be interpreted as the formation and/or expansion of Oxygen Minimum Zones.Finally, in Section 6 we summarise and discuss our results.

The non-spatial system
We consider a conceptual mathematical model that explicitly includes only phytoplankton, zooplankton, and dissolved oxygen.Oxygen is produced by phytoplankton in photosynthesis and consumed by both phytoplankton and zooplankton as needed for their metabolism.In the zero-dimensional (nonspatial) case, the model is given by the following three equations [18,19]: where c, u and v denote, respectively, the concentration of dissolved oxygen, phytoplankton, and zooplankton densities in appropriately chosen dimensionless variables [18,19] and t is dimensionless time.The first term in Eq. (1a) describes the rate of oxygen production in photosynthesis (see [19] for more details) and the second and third terms describe the oxygen consumption by phyto-and zooplankton, respectively, which is assumed to be described by the Monod type kinetics.The first term in Eq. (1b) describes the phytoplankton multiplication; based on earlier work [28,65,68], we consider it as the logistic growth.The second term in Eq. (1b) quantifies the phytoplankton grazing by zooplankton [68], and the third term describes the phytoplankton natural mortality.In Eq. (1c), the first term in the brackets describes the zooplankton growth (with the food assimilation efficacy being assumed to depend on the level of dissolved oxygen [19]), the second term describes the zooplankton natural mortality and the third (quadratic) term describes a combined effect of the competition and the predation by species from higher trophic levels (not included into the model explicitly) [30].Here parameter A is the per capita oxygen production rate, B is the per capita phytoplankton growth rate, σ and µ 1 are natural mortality rates for phyto-and zooplankton, respectively, µ 2 quantifies the strength of nonlinear zooplankton removal.The meaning of other parameters in Eqs.(1) is straightforward.For more details, including biological justification of all terms in the right-hand side of Eqs.(1), see [18,19].
For the reasons mentioned in the introduction, we introduce additional parameter 0 < ε < 1 that quantifies the difference in the phyto-and zooplankton characteristic timescales; in most cases below, we will assume ε ≪ 1.Note that, compared to the original model proposed in [18,19], Eqs.(1) include two essentially new elements, i.e. the quadratic term and a small parameter (cf.Eq. (1c)), which makes model (1)somewhat more realistic.

Steady States Analysis
To explore the dynamics of the temporal model, we study all possible equilibrium points (steady states) of the system (1) and their stability.The system has a total extinction state given by E 0 = (0, 0, 0).To study the dynamics of the system (1) around E 0 we linearize around E 0 and obtain the Jacobian matrix All the eigenvalues of the above matrix J E 0 are real and negative.Therefore, the total extinction state E 0 is always stable.In [19] the authors showed that under some parametric restrictions, the system can have two zooplankton free equilibria.Since the introduction of the quadratic term in (1c) does not affect the zooplankton free equilibrium states, as the zooplankton free state is the form (c, ū, 0) where c is the root of the quartic equation Since the above equation is a fourth-order polynomial, analytical determination of the equilibrium points is nearly impossible.We choose suitable numerical values of the parameters to obtain feasible zooplankton free equilibrium states.Throughout the paper, we fix the parameter values at and suitably vary µ 1 , µ 2 , and ε.For all the values of µ 1 and µ 2 and parameters fixed at (5) the two feasible zooplankton free equilibrium points are given by E 1 = (0.0258, 0.0067, 0) and E 2 = (1.712,2.029, 0).Among the two zooplankton free equilibrium points, one is always a saddle point while the other can be either stable or unstable depending on the parameter values.
For our choice of parameter values ( 5) both E 1 and E 2 are unstable (saddle) with a 2-dimension stable manifold and a 1-dimensional unstable manifold.We plot the c-component of E 1 and E 2 in Fig. 1 with red broken lines.Note the lines are horizontal, because the corresponding steady state values do not depend on µ 1 or µ 2 , as is obvious from Eqs. (3)(4).The system does not possess any other feasible boundary equilibria in this parametric regime.Note that a nontrivial oxygen free equilibrium is not possible in our model (which agrees with biological reasons).Indeed, setting c ≡ 0 in Eq. (1a) immediately leads to u ≡ 0, which in turn leads to z ≡ 0. The coexistence equilibrium of the system (1) is denoted by E * = (c * , u * , v * ) where and c * , u * can be obtained by simultaneously solving the quartic and quadratic equations respectively Because of the complexity of the simultaneous algebraic equations we obtain the coexistence equilibrium points using numerical simulations.The parameter values are fixed at (5) and we choose µ 1 , µ 2 as the bifurcation parameters.Note that for µ 1 > 0, µ 2 = 0, the model was rigorously studied in [19].However, we deal with other cases in this work.For µ 1 = 0, µ 2 > 0, we found that the system (1) has two coexistence equilibrium points which disappears via saddlenode bifurcation for sufficiently smaller values of µ 2 .If we denote the saddle-node bifurcation threshold by µ SL 2 , then for µ 2 > µ SL 2 the system has two coexistence equilibrium points and for µ 2 < µ SL 2 there are no feasible coexistence equilibrium.For µ 1 , µ 2 > 0, the system has a unique coexistence equilibrium throughout the parametric regime of µ 2 , whenever they are feasible.The number of coexistence equilibrium points and their stability is determined by the combination of the parameters µ 1 and µ 2 .The Jacobian matrix evaluated at the coexistence equilibrium E * is and the corresponding characteristic equation is given by The coefficient are where J ij is the (i, j)-th entry of J E * , and ii is the cofactor of J ii .By Routh-Hurwitz conditions the coexistence equilibrium point E * is stable if the following conditions hold p 0 > 0, p 2 > 0 and p 1 p 2 > p 0 .
To observe the change in system's dynamics with the introduction of intraspecific competition among the zooplankton, we choose µ 2 to be the bifurcation parameter.Thus keeping all the parameters fixed, p 1 , p 2 , p 0 are functions of µ 2 only.The coexistence equilibrium E * loses its stability and exhibits oscillatory dynamics through Hopf bifurcation at µ H 2 , which is obtained by solving ).This expression is not mathematically tractable, thus we numerically obtain the Hopf bifurcation threshold µ H 2 = 0.35405 (upto five decimal place) for the parameter values (5), and µ 1 = 0.05, ε = 1.The Hopf bifurcation can be either supercritical or subcritical depending on the magnitude of the parameters µ 1 and µ 2 .For a fixed µ 1 , with an increase in the intraspecific competition (as quantified by µ 2 ), the equilibrium state of the system changes its stability from unstable to stable.The broken line in Fig. 1 represents the unstable branch of the equilibrium and the continuous line represents the stable branch.The red and blue dots on the equilibrium branch are the subcritical and supercritical Hopf bifurcation thresholds, respectively.For a relatively higher value of µ 2 , both the extinction state E 0 and the coexistence state E * is stable, and the system exhibits bi-stability.For µ 1 = 0 the lower branch of the coexistence equilibrium is saddle and it acts as a separatrix between the basin of attraction of E * and E 0 .However, for µ 1 > 0, their basins are separated by the unstable manifold of the saddle boundary equilibrium point.For a fixed µ 2 we observe an increase in the oxygen concentration with an increasing mortality rate of the zooplankton (see Fig. 1).Therefore, an increase in the zooplankton linear mortality rate (as quantified by µ 1 ) and an increase in the nonlinear mortality rate (due to the intraspecific competition and/or predation by higher trophic levels [67], as quantified by µ 2 ) lead to an increase in the stable oxygen level.

Local and global dynamics
The dynamics of the system (1) is determined by its local and global bifurcation structure.However, any comprehensive analytical analysis of the bifurcations is hardly possible for this model due to its algebraic complexity.Instead, we choose µ 1 and µ 2 as bifurcation parameters and obtain two corresponding one parametric bifurcation diagrams; see Fig. 2.These diagrams readily reveal the various Hopf bifurcation scenarios of the coexistence equilibrium.We further show the possibility of the heteroclinic bifurcation and saddle-node bifurcation of limit cycles.In Figure 2, the Hopf bifurcation thresholds (supercritical and subcritical) are marked by black dots and the thresholds for the global bifurcations (heteroclinic and saddle-node bifurcation of limit cycles) are marked by green vertical broken lines.The stable and unstable cycles (and equilibrium) are shown in blue and red colour, respectively.We choose µ 1 or µ 2 as bifurcation parameters to observe the change in the system's dynamics.Note that, the zooplankton free equilibrium points E 1 and E 2 are independent of the bifurcation parameters.Thus, the stability of these equilibrium points does not depend on µ 1 and µ 2 .For our investigation, we choose other parameters as in (5), so that that both E 1 and E 2 are unstable.we plot the c-components of the unstable zooplankton free equilibrium points of the form (c, ū, 0) in a black broken line.
In the absence of intra-specific competition among the zooplankton and neglecting the effect of higher trophical levels (that is, for µ 2 = 0), the system has a unique coexistence equilibrium, which loses its stability through Hopf bifurcation at µ H 1 = 0.398 (considering all the other parameters are fixed at (5)).The Hopf bifurcation is supercritical as the first Lyapunov coefficient is l 1 = −0.358(<0).The unique coexistence equilibrium point is stable for µ 1 > µ H 1 and unstable for µ 1 < µ H 1 .A small amplitude stable limit cycle originates from µ H 1 .The amplitude of the stable cycle increases with decreasing µ 1 till it hits the boundary equilibrium points and disappears via a heteroclinic bifurcation (marked in green).Beyond this threshold, the coexistence equilibrium is unstable and the extinction state E 0 is the only global attractor.On decreasing the mortality rate from µ H 1 the concentration of the zooplankton increases.This increases the oxygen consumption by zooplankton and also decreases phytoplankton production such that beyond a threshold the system collapse.
On the contrary, if we consider the case where the zooplankton mortality is primarily due to the combined effect of the competition and predation by species from higher trophic levels, then one could neglect the natural mortality [67], so that µ 1 = 0. Using other parameter values as in (5), the Hopf bifurcation occurs at µ H 2 = 0.408.Since the first Lyapunov coefficient is l 1 = 1.065(> 0), the Hopf bifurcation is subcritical.Among the two coexistence equilibria as obtained in Fig. 1, one is always saddle, the other is unstable for µ 2 < µ H 2 and stable for µ 2 > µ H 2 .An unstable cycle emerges from the subcritical Hopf bifurcation, increases in size in a narrow domain for µ 2 > µ H 2 , till it disappears via global bifurcation (see Fig. (2)(b)).Beyond this, the stable extinction state E 0 coexists with a stable coexistence equilibrium.We now choose µ 1 = 0.24.Then, the subcritical Hopf bifurcation occurs at µ H 2 = 0.1577.Here the unstable cycle that emerged through Hopf bifurcation is surrounded by a stable limit cycle.The stable and unstable cycle exists in a narrow range of µ 2 and collides at saddle-node bifurcation of limit cycle (µ 2 = 0.1578).On the other hand, for µ 2 < µ H 2 , the unstable equilibrium is surrounded by a stable limit cycle.The amplitude of the stable cycle increases and disappears through heteroclinic bifurcation at µ 2 = 0.15715 (see Fig. (2)(c)).

Slow and fast dynamics
The system (1) for 0 < ε ≪ 1 is referred to as the singularly perturbed system.The timescale parameter ε signifies a clear distinction between two timescales: slow and fast.The change in the zooplankton density occurs at a much slower rate as compared to the change in oxygen and phytoplankton density.We denote the time 't' in system (1) as the fast time and the system (1) is with respect to the fast timescale.
Introducing the slow timescale τ = εt, we obtain the following slow subsystem: We analyze the above slow and fast systems with the help of geometric singular perturbation theory [61,69].The basic idea behind this approach was to decompose the slow and fast systems in its limiting systems (i.e for ε = 0) and study the dynamics of the respective subsystems.In the singular limit ε → 0, we obtain the fast subsystem (oxygen-phytoplankton) of system (1) as follows with v = v 0 (constant zooplankton density).Also, letting ε = 0, in (7) we obtain the slow subsystem or the reduced system as where F, G, H are given above.The set is called critical manifold, which is the collection of all equilibrium points or curves of the fast subsystem.It can be divided into two parts: trivial critical manifold C 0 0 and non-trivial critical manifold The trivial critical manifold is given by the line and the non-trivial critical manifold C 1 0 is the curve of intersection of two surfaces given by The slow subsystem (9) represents the slow change in the zooplankton density over the critical manifold.Considering the slow variable v as the bifurcation parameter, we numerically obtain the critical manifold (black) as shown in Fig. (4).Therefore, along the critical manifold, the zooplankton density changes slowly.Since the two surfaces intersect along a curve, we find the extrema (fold point) of the curve (if any).A point P ∈ C 0 is a fold point of the critical manifold if the fast subsystem exhibits a fold bifurcation.In other words, if we consider the Jacobian matrix of the fast subsystem then at the fold point P , rank(J ) = 1 and the critical manifold is non-hyperbolic.It divides the critical manifold into two halves namely attracting C a 0 and repelling C r 0 sub-manifold.A point p ∈ C a 0 if both the eigenvalues of the matrix J evaluated at p has negative real parts, whereas a point q ∈ C r 0 if atleast one of the eigenvalues of the matrix J evaluated at p has positive real part.To determine the stability of the trivial manifold, we evaluate the matrix J along C 0 0 and thus obtain Since all the parameters involved in the system are positive, the eigenvalues are given by for v ≥ 0. Thus the trivial manifold C 0 0 is stable for v ≥ 0. Fenichel's theorem [69] state that the normally hyperbolic attracting and repelling sub-manifolds, C a 0 and C r 0 respectively, obtained for ε = 0, perturb to locally invariant attracting and repelling sub-manifolds C a ε and C r ε respectively, for ε > 0. Therefore, the dynamics of the full system (1) or ( 7) can be approximated by studying the dynamics of the subsystems obtained for ε = 0.
The slow flow on the critical manifold C 1 0 is given by the slow subsystem (9).We differentiate F (c, u, v) = 0 and G(c, u, v) = 0 implicitly with respect to 'τ ' along the critical manifold, and obtain the dynamics on the critical manifold.This is governed by the following system of equations with suitable initial condition (c 0 , u 0 , v 0 ) ∈ C 1 0 .The slow flow has a singularity whenever G c F u − G u F c = 0, which holds at the fold point P. Thus, the solution blows up at this point.Whenever the fold point is called the jump point, and the trajectory jumps from the proximity of the fold point to another attracting critical manifold.However, when both the fold point is called canard point.At this point, the trajectory can pass through the proximity of the fold point and follow the repelling manifold for O(1) time.The coexisting equilibrium loses its stability and exhibits oscillatory dynamics through Hopf bifurcation (supercritical or subcritical), which we discussed in the previous section.In a classical slow-fast setting, the small cycles originating from Hopf bifurcation bifurcates to canard cycles (with or without head) and further to relaxation oscillation, thus exhibiting canard explosion [70].These cycles are composed of slow and fast segments, where the slow flow occurs along both attracting and repelling sub-manifold of the critical manifold.It exhibits fast flow when the trajectory jumps to another stable portion of the critical manifold.However, for the system (1), small canard cycles (without head) emerge from Hopf bifurcation.The amplitude of the cycles increases in a narrow parametric range, eventually leading to complete extinction.We prove this fact in the following theorem.
Theorem 4.1.Assume the fold point P is a canard point for ε > 0. Then the system (1) has small amplitude canard cycles (without head) originating from Hopf bifurcation but there does not exist any relaxation oscillation.
Proof.In the slow-fast setting, we denote the Hopf bifurcation as singular Hopf bifurcation since the eigenvalues of the Jacobian matrix of the system (1) evaluated at E * has purely imaginary complex eigenvalue of the form λ 1,2 = ± iω(µ H 2 , ε) such that lim ε→0 ω(µ H 2 , ε) = 0.The singular Hopf bifurcation occurs at O(ε) from the fold point P. We assume that the fold point is the canard point.The small limit cycle originating from Hopf bifurcation grows in size through a sequence of canard cycles.With the decrease in ε, the amplitude of the cycle increases, and after a certain threshold, the trajectory jumps from the vicinity of the fold point P close to C 0 0 .The equilibrium point E 0 lies on the trivial critical manifold.The eigenvalues of the Jacobian matrix J E 0 are −1, −σ, −εµ 1 and the corresponding eigenvectors are (1, 0, 0) , 1, 1 − σ A , 0 , and (0, 0, 1) .
Therefore the critical manifold C 0 0 coincides with the eigenvector.Thus, any trajectory on C 0 0 converges to E 0 .We cannot construct any singular orbit consisting of concatenated slow segments on C 1 0 and C 0 0 , and fast fibers while leaving the respective manifolds.Hence, the global return mechanism, which is necessary for the existence of classical relaxation oscillation, fails as all the trajectories converge to the stable equilibrium E 0 .The horizontal red (broken) line shows the steady state value of the coexistence equilibrium when it is unstable.The vertical black (broken) line marks the singular Hopf bifurcation threshold (occurring at µ 2 = 0.1007) and the vertical green (broken) line at µ 2 = 0.099171 indicate the threshold for the system collapse (plankton extinction and oxygen depletion).Other parameters of the system are given in (5).
We illustrate this phenomenon with the help of a numerical example in Fig. 3.We use the parameter values as in (5) along with ε = 0.5.Consider a hypothetical value µ 1 = 0.24.The effect of ε on the dynamics of the system can be observed if we compare the Fig. 3(a) with Fig. 2(c).The singular Hopf bifurcation occurs at µ H 2 = 0.1638 is subcritical (l 1 = 0.037).Small unstable canard cycles emerge from the canard point P (1.24, 1.01, 0.89), which grows in size with a slight increase in µ 2 (Fig. 3 (a)).A large amplitude stable cycle emerges from a heteroclinic bifurcation that coexists with the stable equilibrium, separated by an unstable canard cycle.We observe that the size of the stable cycle shrinks in an extremely narrow parameter interval and disappears at a saddle node bifurcation of limit cycles.
We now consider µ 1 = 0.3 and ε = 0.5, then the singular Hopf bifurcation occurs at µ H 2 = 0.1007.The first Lyapunov coefficient is l 1 = −1.1818,hence the singular Hopf is supercritical.From the canard point P (1.255, 1.035, 0.898), small stable canard cycles originate.We show the change in the amplitude of the canard cycles with decreasing values of µ 2 in Fig. 3(b).This depicts that the system becomes unstable with a decrease in the strength of the intra-specific competition among the zooplankton.The transition from the stable, steady state to the oxygenfree state, indicating complete population collapse, takes place in an extremely narrow interval of the rate of intraspecific competition.That is, for µ 2 ∈ (0.099171, 0.1007).At µ 2 = 0.09917, when the size of the limit cycle explodes, the trajectory converges to the origin, which is illustrated in Fig. 4(a).The time series of the trajectory is shown in Fig. 4(b).This implies that the system cannot further sustain any large amplitude oscillations.The rise in the amplitude of the phytoplankton level beyond a threshold can act as an indicator of population collapse.The system can therefore be driven to total extinction by pushing it far enough to reach the fold point.Along with µ 2 , the timescale separation plays a critical role in identifying this narrow parametric regime.For a larger timescale separation (i.e.smaller ε), the coexistence of oxygen-plankton occurs in a significantly narrower interval, and the system is more vulnerable to perturbation.

The spatial system
In the real world ocean environment, the spatial distribution of both plankton and dissolved oxygen is remarkably heterogeneous, sometimes showing the variability by an order of magnitude of even more [27,37,50,51,71].Correspondingly, in this section, we consider a spatially extended model (1) where the oxygen concentration and the phytoplankton and zooplankton densities vary with both time and space.The ocean system is three-dimensional; however, in this paper, for the sake of simplicity, we only consider one horizontal spatial dimension.We regard it as the position along the ocean surface.In terms of ocean observations, it corresponds to a transect across the study area.
Trying to keep the model as simple as possible, we avoid using explicit dependence on the vertical dimension (i.e. the depth).Correspondingly, we use the well-mixed layer approximation [28,65,68] to assume that the vertical distribution of plankton and oxygen is approximately uniform within the photic (upper) ocean layer where most of photosynthetic oxygen production takes place.
The transport of any substance in the ocean takes place primarily due to the water movement.The movement in the horizontal direction occurs either due to an ocean current or marine turbulence (or their combination).Here we chose to focus on the effect of unbiased (isotropic) movement, hence taking into account only the effect of turbulence, which we describe as the turbulent diffusion quantified by a certain diffusion coefficient [72,73].
We, therefore, arrive at the following equations: Here c(x, t), u(x, t), and v(x, t) are, respectively, the oxygen concentration and the phytoplankton and zooplankton densities at the (horizontal) location x and time t; D c , D u and D v denote the diffusion coefficients for oxygen, phytoplankton, and zooplankton.Note that phytoplankton and the dissolved oxygen can be regarded as a 'passive substance', i.e. their spatial transport is entirely determined by the water flows; hence D c = D u = D T where D T is the turbulent diffusion coefficient.However, zooplankton has a certain ability to self-motion.Combined with the effect of turbulent mixing, it can result in a value of D v ̸ = D T .Whether D v is larger or smaller, then depends on the zooplankton movement pattern.In case zooplankton movement is entirely random (e.g. can be regarded as Brownian motion), then one can expect that D v > D T .
In case zooplankton exhibits a homing behavior, then it is likely that D v < D T .
Along with the temporal parameters, we now non-dimensionalize the space as x = x √ Dc .Removing the tilde for the simplicity of the notation, we obtain the following dimensionless spatial model: where D = Dv Dc .Equations ( 14) are considered inside the spatial domain Ω = {x ∈ (0, L)} where L is thus the length of the domain.
Equations ( 14) must be complemented with the initial conditions, which we consider in the following form: That is, at t = 0 the steady state densities are perturbed within a small area at the center of the domain Ω.
For the boundary conditions, we consider the zero-flux conditions: The model is solved numerically with L = 500.We use the Euler method for the temporal part and five points central difference scheme for the diffusion part along with ∆x = 1 and ∆t = 0.01.

Turing instability
To study the spatial distribution of oxygen and plankton, we start our analysis in a neighborhood of a homogeneous steady-state solution of ( 14).Time-independent or a steady state solution (c(x), u(x), v(x)) of the system ( 14)-( 16) satisfies the following system of equations The coexistence equilibrium E * of system (1) corresponds to a homogeneous steady state solution of the above system.To study the dynamics of the above system near the homogeneous steady-state solution E * , we give small heterogenous perturbation as with 0 < ξ 1 , ξ 2 , ξ 3 ≪ 1.The parameter k is the wavenumber of the eigenfunction, and λ is the eigenvalue determining the temporal growth of the corresponding k th mode.We obtain the linearized system as where Z ≡ (z 1 , z 2 , z 3 ), D ≡ diag(1, 1, D).For the non-trivial solution of the above system (19), the eigenvalues λ are determined by the roots of the characteristic polynomial det(λI − J E * + Dk 2 ) = 0, which is written explicitly as where ) and J ij and ii are the same as obtained during the analysis of the temporal part.We, therefore, obtain the necessary and sufficient conditions for Turing instability as p 2 (0) > 0, p 0 (0) > 0, p 1 (0)p 2 (0) > p 0 (0) and p 0 (k 2 ) < 0, for some k. ( Therefore, the Turing instability occurs at a critical wave number, k = k T , where p 0 (k 2 ) achieves a minimum and p 0 (k 2 T ) = 0.This gives where where J ij are the elements of the Jacobian matrix, cf.Eq. ( 6).Because of the complexity of the expression and the large number of parameters involved, the critical wave number k T corresponding to the Turing bifurcation has to be computed numerically.For a feasible k T , the model describes the formation of spatial patterns, as is shown below.

Impact of diffusion on Oxygen Minimum Zone
We now look into the spatio-temporal dynamics of the system ( 14) with the initial conditions (15).Our goal is to reveal typical dynamical regimes (in particular, pattern formation scenarios, if any) for parameters inside and outside of the Turing domain.In a pure Turing domain all the Turing instability conditions as discussed above holds.However, in a Turing-Hopf domain, the homogeneous steady state is unstable under both temporal and spatio-temporal perturbations.
We fix the parameter values as in ( 5) and let µ 1 = 0, µ 2 = 0.41 and consider different values of the diffusivity ratio D. Fig. 5 shows typical patterns in the distribution of oxygen for different diffusivity rates.Here Fig. 5a,b and Fig. 5c are obtained, respectively, for parameters outside and inside the Turing instability domain.We notice that, in all three cases, the evolution of the initial conditions soon leads to the formation of a patch where the oxygen concentration is much lower than its steady state value.We interpret this dynamics as the formation of an OMZ.Further evolution of the emerging OMZ can be significantly different depending on D. When the diffusivity ratio is small, i.e.D < 1, the OMZ created at the early stage grows with time and eventually spreads over the entire domain.Interestingly, the growing OMZ has a fine structure.A closer look at the dynamics shown in Fig. 5a reveals that, at any time t during the transient stage of the OMZ expansion, it consists of three or four subdomains with very low oxygen level separated by narrow spatial intervals where the oxygen level is larger than its steady state value c * .This fine structure disappears after the expanding OMZ hits the domain boundaries; at a later time the oxygen level is low over the entire domain, which can be interpreted as the global anoxia.
An increase in the diffusivity ratio to D = 1 results in a qualitative change in the dynamics; see Fig. 5(b).In this case, the OMZ formed at an early stage of the system dynamics show almost no spatial growth remaining localised around the centre of the domain.At any spatial position inside the OMZ, the oxygen level distinctly oscillates with time between a very low level (approximately 0.1c * ) to a high level (of about 1.3c * ).
A further increase in D leads to another qualitative change in the dynamical pattern.Figure 5(c) shows the results obtained for D = 5.In this case, the system satisfies the Turing instability condition with the critical wavenumber k 2 T = 0.1095.The evolution of the initial conditions leads to the formation, inside a certain subdomain, of a periodic spatial pattern where low oxygen patches alternate with high oxygen patches.The subdomain containing this periodic structure grows with time and eventually occupies the whole domain, so that at a large time the periodic spatial distribution becomes stationary.
We note here that the pattern of the OMZ formation and spread is relatively robust to the initial conditions.For instance, if at t = 0 the spatial distribution of oxygen is perturbed along with that of phytoplankton, the emerging patterns are similar to the ones shown in Fig. 5.One example is shown in Fig. 6.In this case, the initial conditions (15) are slightly modified, so that, at the center of the domain both c(x, 0) and u(x, 0) are less than their steady state values.
It is readily seen that the top and the bottom of Fig. 6 show very similar patterns, with the only difference that the spatial size of the emerging OMZ becomes larger for the modified initial conditions.

Impact of timescale separation on Oxygen Minimum Zone
From the mathematical analysis in the subsection 5.1, we obtain that the critical wavenumber k T for Turing instability depends on the timescale separation ε, as the elements J 31 , J 32 and J 33 of the Jacobian matrix in Eq. 23 depend on ε.Thus, one can expect that the boundaries of the parameter ranges where the Turing and the Turing-Hopf instability occur (resulting in pattern formation) can shift with a decrease in ε, i.e. with an increase in the timescale separation.Numerical simulations confirm that this is indeed the case.Having fixed µ 1 = 0, µ 2 = 0.41, D = 5 and other parameters as in ( 5) and only varying ε, we obtain stationary patterns in both the Turing domain and in the Turing-Hopf domain for ε ≥ 0.18.A typical pattern is shown in Fig. 5(c).With a decrease in ε, the emerging stationary pattern has a similar nature of alternating patches of high and low oxygen level as for ε = 1 (cf.Fig. 5c) but the size of the patches becomes larger; e.g.see Fig. 7(a).Also, the emergence of the stationary periodic pattern is preceded by rather long transient dynamics when the oxygen concentration and the plankton densities exhibit irregular oscillations (see Fig. 7(d)).
With a further decrease in ε, the dynamics becomes qualitatively different.The emerging pattern is not spatially periodic any more; see Fig. 7(b).Apart from the large OMZ formed around the centre of the domain at the early stage of system's dynamics, there are two large OMZs at the sides of the domain.At a later time, these patches of low oxygen level break to a number of smaller patches of variable size.The dynamics is not becoming stationary at any time as the patches keep changing their size (and some of them also their location).The dependence of spatially average densities is distinctly irregular (see Fig. 7(e)) suggesting chaotic dynamics.This kind of dynamic pattern is observed for 0.07 < ε < 0.18.
With a further decrease in ε (below ε = 0.07), the system's dynamics undergo a regime shift.For ε < 0.07, the transient apparently chaotic dynamics only last for a finite time.After a sufficiently long time, the system experiences a catastrophic change when over a short transition   time the oxygen concentration fast drops to a very small value (and eventually to zero) over the entire spatial domain.An example of such regime shift is shown in Fig. 7(c,f) (obtained for ε = 0.06).The entire area becomes a dead zone (with low or no oxygen), which can be interpreted as the global anoxia.Along with oxygen, the phyto-and zooplankton densities go to zero as well (after several irregular oscillations of increasing amplitude, cf.Fig. 7(f)), obviously signifying their extinction.For both µ 1 > 0 and µ 2 > 0, the system's dynamics becomes different and exhibits a somewhat greater variety of dynamical regimes.As one example, Figs. 8(a,e) shows the spatiotemporal dynamics for µ 1 = 0.24, µ 2 = 0.1575, D = ε = 1 and other parameters the same as in Fig. 5(b,e).It is readily seen that, in this particular case, the evolution of the initial condition does not lead to formation of OMZ.It only leads to small fluctuations in the oxygen level and plankton densities around the location of the initial perturbation, with the spatial distribution being uniform in the rest of the domain.
A decrease in ε leads to the emergence of the OMZ.It first appears at the position of the initial perturbation (i.e.near the centre of the domain); see the bottom of Figs.8(b,c).At a later time, it breaks to several patches that fast spread over the entire domain.The spatiotemporal dynamics is apparently chaotic for ε = 0.5 but becomes more regular for ε = 0.25, cf.Figs.8(f,g).
A further decrease in ε below a certain critical value results in a regime shift, e.g.see Figs. 8(d,h) obtained for ε = 0.2.In this case, a large OMZ is formed at an early stage of system's dynamics (see the bottom of Fig. 8(d)).However, after a relatively short time the oxygen concentration fast drops to zero over the entire domain: the global anoxia occurs accompanied by the plankton extinction.

Conclusion
Over the last few decades, there has been growing evidence of a decline of the dissolved oxygen concentration in the ocean [6,7].This has not only been recognized as a catastrophic threat to the marine ecosystems [55] but also as a potential threat to mankind [49] and to terrestrial ecosystems, as marine phytoplankton contributes about 70% to the total atmospheric oxygen.Any significant decline in the global phytoplankton abundance and/or a decrease in the oxygen production rate in phytoplankton photosynthesis will inevitably lead to a decline in the global stock of the atmospheric oxygen [18,47].Thus, marine ecosystems, phytoplankton in particular, play a crucial role in maintaining the habitable Earth [23].
In spite of the apparent importance of the above issues, mathematical models addressing the change in the oxygen concentration as a component of the coupled phytoplankton-oxygen dynamics are rare in the literature.As one exception, a generic three-component oxygen-phytozooplankton model was developed in [19] (and further investigated in [18,20,74]).It has been shown that the formation of areas with a low oxygen concentration (i.e.OMZs) is in fact an inherent property of the self-organised plankton-oxygen spatiotemporal dynamics, but it can be exacerbated by the effect of global warming, potentially leading to global anoxia.
The model developed in [18,19], however, missed several important features of the marine ecosystem's dynamics, hence making the prediction of emerging global anoxia somewhat questionable.In this paper, we have considered a nontrivial extension of the original model that includes factors such as zooplankton inherent competition, cannibalism and/or the effect of zooplankton's consumers from upper trophic layers, e.g.fish.Another important factor is the existence of different timescales for phyto-and zooplankton growth, as the latter is usually much slower than the former.
The properties of the extended model have been analysed in much detail using a combination of analytical and numerical tools.We first consider a non-spatial version of the model described by a system of three nonlinear ODEs (for oxygen, phytoplankton, and zooplankton, respectively) to study the variation in oxygen level and plankton densities over time.In addition to the results earlier obtained in [19,46], we have shown that the number of coexisting steady states depends both on the linear mortality rate (µ 1 ) and the rate of zooplankton intraspecific competition/consumption (quantified by coefficient µ 2 ).For µ 1 = 0, there exists two feasible coexisting steady states where the lower oxygen level is always unstable, and the higher state changes its stability with increasing µ 2 .Whereas for µ 1 ̸ = 0, we obtain a unique feasible steady state which changes its stability from unstable to stable for a strong intraspecific competition.Along with this, the extinction state is always stable.We also found that an increase in the rates of zooplankton linear mortality (µ 1 ) and nonlinear mortality (µ 2 ) leads to an increase in the oxygen abundance.
In order to better understand the relative importance of the linear and nonlinear mortality and their effect on the temporal dynamics of the system, we have considered three cases: (a) µ 1 ̸ = 0, µ 2 = 0, (b) µ 1 = 0, µ 2 ̸ = 0, and (c) µ 1 ̸ = 0, µ 2 ̸ = 0.Because of the complexity of the system, this has mostly been done through numerical simulations.For case (a), the unique steady state is stable for higher values of µ 1 , and it loses its stability through supercritical Hopf bifurcation.Small stable cycles originate, increasing its size in a small interval of µ 2 .Beyond that, the system cannot further withstand an increase in the amplitude of the cycle leading to complete collapse (see Fig. 2a).However, for case (b), the Hopf bifurcation is subcritical, and the system converges to stable steady state for higher values of µ 2 .In this case, an unstable cycle is formed (see Fig. 2b).Case (c) is a combination of the above two cases.Here, in a very narrow domain, the system exhibits tri-stability, with a stable steady state, a stable cycle, and the extinction state.The two cycles appear through saddle-node bifurcation of limit cycle, and the disappearance of the unstable cycle is through subcritical Hopf, and that of the outer stable cycle is through heteroclinic bifurcation (see Fig. 2c).
Having analysed the effect of different timescales (cf."slow-fast system"), we obtained the critical manifold of the slow subsystem.The extremum (the fold point) of the critical manifold acts as an extinction threshold: if the system is pushed beyond the fold point (e.g. by the choice of the initial conditions), the dynamics will eventually lead to plankton extinction and oxygen depletion, although the extinction/anoxia can be preceded by a long period of oscillations (cf.Fig. 4).Note that a decrease in the nonlinear mortality rate µ 2 has a similar effect on the system's persistence.A decrease in µ 2 below a certain critical value first destabilises the coexistence steady state resulting in oscillatory dynamics (see Fig. 3).A further decrease (below another critical value) leads to the canard explosion.However, no relaxation oscillations emerge (see Theorem 4.1); instead, the system's trajectory goes to the origin, which obviously corresponds to the extinctions and anoxia.
We then considered the spatially explicit system to study the distribution and spatiotemporal dynamics of oxygen and plankton.The spatially explicit model consists of three reactiondiffusion equations where the diffusion terms account for the effect of lateral turbulent mixing for the dissolved oxygen and phytoplankton and for the combined effect of turbulence and self-motion for zooplankton.Note that, because of the latter, zooplankton diffusivity can be expected to differ from that of phytoplankton and oxygen.Moreover, it can differ significantly.The interplay between the ordinary Fickian diffusion (in our case "biodiffusion" resulting from zooplankton random movement) and the turbulent mixing is known to be highly nonlinear [73,75].The ordinary diffusion, although itself often being orders of magnitude less intensive than the turbulent mixing, accelerates the turbulent diffusion significantly [73,75].In turn, for the diffusivity ratio being greater than one, the model can exhibit pattern formation due to the Turing instability; see Sections 5.1 and 5.2.
Due to its mathematical complexity, the spatially explicit reaction-diffusion model is not analytically tractable.Therefore, we have investigated its properties through extensive numerical simulations, with a special attention to regimes that result in the formation of patterns containing areas with low oxygen level and/or regimes resulting in global oxygen depletion.Using the initial condition as a localised perturbation of the spatially uniform steady state, we have obtained that the system dynamics typically lead to the formation of strongly heterogeneous spatial distribution that includes one or several areas (patches) with a very low oxygen level, which we interpret as the formation of OMZ.Interestingly, the patterns emerge both in and outside of the Turing domain and hence, for different parameter values (e.g. the diffusivity ratio being larger or smaller than one) can be attributed to different dynamical mechanisms, i.e.Turing or non-Turing.Except for some rare cases (cf.Fig. 5b), the OMZ formed at an early stage of the system dynamics fast spreads over the entire domain, often generating multiple patches, e.g.see Figs. 5(a,c), 7(a,b,c) and 8(b,c,d), the size and number of the emerging smaller OMZs varying with the parameter values.
The spread of the emerging pattern (a mixture of patches with high and low oxygen level) can lead to a different outcome.It can result in a self-sustained pattern, which, in the large time limit, can be stationary (cf.Figs.5c and 7a) or dynamic (Fig. 7b).Alternatively, it may eventually lead to an unsustainable pattern -a regime shift -when, after a certain time, the oxygen concentration fast drops to very small values over the entire domain (cf.Fig. 5a, 7c and  8d).Arguably, this may be interpreted as the onset of the global anoxia.
Note that there is a subtle interplay between the zooplankton linear mortality rate µ 1 and the difference in the timescales.In the special case of the same timescales (ε = 1), the effect of a non-zero zooplankton linear mortality makes the system somewhat more sustainable: while an initial perturbation leads to the formation of deoxygenated patch (OMZ) at the center of the domain in case µ 1 = 0 (Fig. 5b), it only leads to small fluctuations in the oxygen level in case µ 1 > 0 (Fig. 8a).However, for ε < 1 and with a further increase in the timescale separation (i.e. for a smaller ε), the effect of zooplankton mortality becomes rather opposite making the system less sustainable.For instance, for µ 1 > 0 the global anoxia occurs already for ε = 0.2 (see Fig. 8d) but for µ 1 = 0 the dynamics remains sustainable (i.e.no global anoxia) for ε = 0.1.
For a fixed value of µ 1 , an increase in the timescale separation alone makes the dynamics less sustainable.A decrease in ε tends to lead to a larger size of the initially formed OMZ; e.g.see Fig. 8(b,c,d), eventually resulting in global anoxia and extinctions when ε becomes sufficiently small.That happens both for µ 1 = 0 and µ 1 > 0, cf.Figs.7(c) and 8(d), although the succession of spatial patterns preceding the onset of anoxia is different between the two cases.Apparently, our study leaves open questions.Firstly, recall that our model is conceptual; it only takes into account the interaction between oxygen and plankton but not with other components of the complicated marine food web.It has been shown in [18] that, in case of a trophic chain, the effect of higher trophic levels only makes the regime shift -the catastrophe of oxygen depletion and plankton extinction -more likely as the three-component model (1) provides an upper bound for a longer trophic chain.An open question however remains as to how the dynamics may change in case of a web rather than chain, for instance to account for effect of bacteria or detritus.Secondly, the description of turbulent mixing as the turbulent diffusion is somewhat simplistic; in particular, it completely disregards the fact that the turbulent mixing is multiscale and nonlocal [73,75].Although the model with the turbulent diffusion is arguably a sensible first step, a more advanced approach should involve a more realistic description of turbulence.These issues will become a focus of future work.

Figure 1 :
Figure 1: The plot of c-component of the coexistence equilibrium for the parameter values given in (5) as a function of parameter µ 2 for a few different values of µ 1 : µ 1 = 0 (blue), µ 1 = 0.05 (magenta), µ 1 = 0.25 (green), and µ 1 = 0.3 (black).The black dot represents the saddle-node bifurcation threshold (c SL , µ SL 2 ).The inset shows the nature of the equilibrium branch near the saddle-node bifurcation threshold in log-log scale.The red horizontal lines represent the c-component of the two zooplankton-free equilibrium points.The red dots on the equilibrium branches represent subcritical Hopf bifurcation threshold, while the blue dots on the equilibrium branches represents supercritical Hopf bifurcation threshold.

Figure 2 :
Figure 2: Bifurcation diagram of the system (1) at the parameter values (5) for different combination of µ 1 and µ 2 .The stable equilibrium and limit cycles are marked in blue.The unstable cycle and coexistence equilibrium are marked in red (continuous and broken respectively).The c-component of the boundary equilibrium points is shown by broken black lines.The global bifurcation thresholds are marked in vertical broken green line.The black dot represents the Hopf bifurcation threshold (subcritical and supercritical).

Figure 3 :
Figure 3: The change in the amplitude of the canard cycle emerging from the singular Hopf bifurcation with varying µ 2 is shown for (a) µ 1 = 0.24, ε = 0.5, and (b) µ 1 = 0.3, ε = 0.5.The blue lines show the steady state value of of the coexistence equilibrium when it is stable and the maximum and minimum amplitude of the stable canard cycle when the equilibrium is unstable.The horizontal red (broken) line shows the steady state value of the coexistence equilibrium when it is unstable.The vertical black (broken) line marks the singular Hopf bifurcation threshold (occurring at µ 2 = 0.1007) and the vertical green (broken) line at µ 2 = 0.099171 indicate the threshold for the system collapse (plankton extinction and oxygen depletion).Other parameters of the system are given in(5).

Figure 4 :
Figure 4: (a) A trajectory (blue) converging to the origin (extinction) after oscillations of increasing amplitude obtained for µ 1 = 0.3, µ 2 = 0.09917 and ε = 0.5.The other system parameters are mentioned in equation (5).The surfaces F = 0 and G = 0 are shaded in green and brown, respectively.The black curve on the intersection of these surfaces is the critical manifold C 1 0 .The single and double arrows represent slow and fast motion, respectively.(b) The corresponding dependence of the phytoplankton density on time.

Figure 5 :
Figure 5: Transition of spatio-temporal dynamics of oxygen from non-Turing (panels (a) and (b)) to Turing (panel (c)) pattern formation for µ 1 = 0, µ 2 = 0.41, ε = 1 and different values of D. All other parameters are given in (5); the corresponding steady state value c * ≈ 1.2.The auxiliary red lines help to reveal the properties the oxygen distribution of oxygen at a given moment of time (as in panels (a) and (c)) or at a given location in space (as in panel (b)); see details in the text.

Figure 6 :
Figure 6: (Top) Zoomed plot of Fig.5(a), which is obtained using the initial condition (15).(Bottom) The initial conditions are of the form (15) but with c(x, 0) = c * − 0.5 and u(x, 0) = u * − 0.2 for |x − L 2 | < 10.All the parameters are same as in Fig. 5(a).The modified initial conditions therefore result in the formation of the OMZ of a larger size.