Pattern Formation in a Model Oxygen-Plankton System

Decreasing level of dissolved oxygen has recently been reported as a growing ecological problem in seas and oceans around the world. Concentration of oxygen is an important indicator of the marine ecosystem’s health as lack of oxygen (anoxia) can lead to mass mortality of marine fauna. The oxygen decrease is thought to be a result of global warming as warmer water can contain less oxygen. Actual reasons for the observed oxygen decay remain controversial though. Recently, it has been shown that it may as well result from a disruption of phytoplankton photosynthesis. In this paper, we further explore this idea by considering the model of coupled plankton-oxygen dynamics in two spatial dimensions. By means of extensive numerical simulations performed for different initial conditions and in a broad range of parameter values, we show that the system’s dynamics normally lead to the formation of a rich variety of patterns. We reveal how these patterns evolve when the system approaches the tipping point, i.e., the boundary of the safe parameter range beyond which the depletion of oxygen is the only possibility. In particular, we show that close to the tipping point the spatial distribution of the dissolved oxygen tends to become more regular; arguably, this can be considered as an early warning of the approaching catastrophe.


Introduction
Oxygen is an important indicator of the marine ecosystem's health. Its depletion may lead to severe ecological problems resulting in mass mortality of marine fauna. An understanding of the dynamics of the dissolved oxygen is therefore important in order to forecast emerging anoxic events and the corresponding ecological disasters. For this reason, the dynamics of oxygen in marine ecosystems has been a focus of both empirical and theoretical research for a few decades [1][2][3][4][5]. In particular, Marchettini et al. [4] suggested a 'realistic' multi-component model of biochemical processes of a lagoon system that explicitly includes oxygen. Using a simpler 'conceptual' model, Allegretto et al. [1] addressed the possibility of periodic fluctuations in the oxygen concentration considering an Italian coastal lagoon as a paradigmatic real-world marine ecosystem. An oxygen-algae model was introduced in Misra et al. [5] to study the oxygen depletion under the effect of some controlling external factors. In another work by the same group, the effect of the time delay due to detritus re-mineralization on oxygen depletion was studied [6]. Hull et al. [2] investigated dissolved oxygen concentration in another multi-component system paying particular attention to the role of bacteria and the effect of the environmental forcing (through wind, solar radiation and temperature) on lagoon dynamics. In yet another study, the seasonal and daily dynamics of the dissolved oxygen in Mediterranean coastal lagoons was investigated [3]. Caruso et al. [7] considered a stochastic model linking the plankton dynamics to the climatic changes and the variations in 18 O isotope.

Mathematical Model
The structure of the model is shown schematically in Figure 1. Flows of matter through the system are indicated by arrows, and the quantification of these flows are given on each arrow respectively. This block-chain model therefore describes the increase in oxygen concentration as a result of photosynthetic activity by phytoplankton. The produced oxygen is consumed in metabolic activities involving respiration for both phytoplankton and zooplankton. Furthermore, phytoplankton is grazed upon by zooplankton. Therefore, the role of zooplankton is twofold: to control phytoplankton density by predation and to consume oxygen through breathing. Most of the theoretical studies, however, considered the oxygen dynamics using the approximation of the "well-mixed layer", hence almost completely disregarding the effect of spatial structure and, correspondingly, using non-spatial or spatially-implicit models. Meanwhile, it was observed in several empirical studies that the spatial distribution of dissolved oxygen is strongly heterogeneous often resulting in a prominent patchy structure [8,9]. The effect of the spatial heterogeneity on the oxygen dynamics remains poorly investigated. In our recent work [10][11][12], we considered a novel model of coupled plankton-oxygen dynamics and showed that it is capable of producing transient, irregular, self-organized spatial patterns. It was shown that the pattern formation can affect the limits of the system's stability with respect to environmental forcing such as the global warming [10]. Moreover, it was also shown that the properties of the spatiotemporal pattern evolve when the system is approaching its stability limits, so that the change in the pattern can be regarded as an early warning signal of the approaching tipping point [10].
The above studies, albeit providing an insight into the effect of space on the oxygen dynamics in a marine ecosystem, only considered pattern formation in the plankton-oxygen system with one spatial dimension. Since the properties of the environment were assumed to be uniform (the model parameters being constant, i.e., not depending on space), arguably the results can be interpreted as a transect along the ocean surface. Those studies therefore left open the question what the dynamics can be in a more realistic 2D system. The goal of this paper is to address this issue. We consider a spatially explicit model of the coupled plankton-oxygen dynamics in two spatial dimensions. Mathematically, the model is given by a system of partial differential equations of the reaction-diffusion type. We show that it exhibits a rich spatiotemporal dynamics resulting in the formation of a variety of patterns. The environment is assumed to be uniform, hence the emerging patterns are self-organized.

Mathematical Model
The structure of the model is shown schematically in Figure 1. Flows of matter through the system are indicated by arrows, and the quantification of these flows are given on each arrow respectively. This block-chain model therefore describes the increase in oxygen concentration as a result of photosynthetic activity by phytoplankton. The produced oxygen is consumed in metabolic activities involving respiration for both phytoplankton and zooplankton. Furthermore, phytoplankton is grazed upon by zooplankton. Therefore, the role of zooplankton is twofold: to control phytoplankton density by predation and to consume oxygen through breathing. The structure of the conceptual model describing the interactions between oxygen (c), phytoplankton (u) and zooplankton (v). Arrows show the flows of matter through the system, and the parameterizations of the rates are as labelled. Phytoplankton produces oxygen through photosynthesis during the day-time then consumes it during the night [13]. Zooplankton feeds on phytoplankton and consumes oxygen through breathing; see further details in the text. The structure of the conceptual model describing the interactions between oxygen (c), phytoplankton (u) and zooplankton (v). Arrows show the flows of matter through the system, and the parameterizations of the rates are as labelled. Phytoplankton produces oxygen through photosynthesis during the day-time then consumes it during the night [13]. Zooplankton feeds on phytoplankton and consumes oxygen through breathing; see further details in the text. We begin with the nonspatial system (for the moment still using the approximation of a well-mixed upper layer) which is described by the following equations [11,12]: where c is the concentration of the dissolved oxygen, and u and v are the densities of phytoplankton and zooplankton, respectively. The term A f (c) describes the rate of oxygen production per unit phytoplankton mass where f (c) describes the rate of increase in the concentration of the dissolved oxygen due to its transport from phytoplankton cells to the surrounding water and parameter A quantifies the rate of oxygen production inside the cell; A is known to depend on environmental factors such as the availability of sunlight, the availability of CO 2 and the water temperature. The terms u r and v r describe, respectively, the rates of the oxygen consumption by the phytoplankton (as required for its metabolism) and zooplankton (as needed for its breathing). The last term in Equation (1) accounts for the losses of oxygen due to its natural depletion (e.g., due to biochemical reactions in the water) with rate m. Function g(c, u) is the per capita phytoplankton growth rate. We consider it to be a function of the oxygen concentration as, apparently, a certain minimum amount of oxygen is needed for the phytoplankton metabolism, and there is biological evidence for that [14]. Function e(u, v) describes the grazing of zooplankton on phytoplankton, i.e., predation. The consumed phytoplankton biomass is transformed into the zooplankton biomass with efficiency ϕ; see the first term in the right-hand side of Equation (3). Since the well-being of zooplankton obviously depends on the oxygen concentration (so that, ultimately, it dies if there is not enough oxygen to breathe), we assume that ϕ = ϕ(c). The terms σu and µv describe the natural mortality of phytoplankton and zooplankton with rates σ and µ, respectively. Using some available data along with relevant biological argument, all functional responses in Equations (1)-(3) can be specified; the details of how such a model is derived are described in detail in a previous publication [12]. Having then introduced dimensionless variables [10,12], we arrive at the following system: where B and β are, respectively, the phytoplankton linear growth rate and the efficiency of the mass turnover in the large oxygen limit (i.e., when oxygen is not a limiting factor), γ is the competition rate (e.g., due to the self-shading when the plankton density is large), and h, h 1 , h 2 , h 3 , h 4 are half-saturation constants. All coefficients are dimensionless. Note that Equations (5) and (6), if considered separately from Equation (4) (hence neglecting oxygen dynamics altogether and keeping the oxygen concentration constant at some hypothetical steady value), coincide with the classical Rosenzweig-McArthur predator-prey model. However, we mention here that a formal reduction of the System (4)-(6) to the two-component phyto-zooplankton model is not possible as dc/dt = 0 cannot give a steady state value of c in case of changing variables u and v. System (4)-(6) has a few steady states [12]: the extinction state (0, 0, 0), the coexistence state (c 0 , u 0 , v 0 ), and two boundary 'zooplankton-free' states (c 1 , u 1 , 0) and (c 2 , u 2 , 0). Whilst the extinction state exists for any parameter value, the nontrivial steady states only exist in a certain parameter range; see [10,12] for details.
Equations (4)-(6) do not contain space. However, in case the water layer is not well mixed, the plankton density and the oxygen concentration become functions of space. In order to make the model spatially explicit, we need to include the transport processes that are responsible for carrying around the system's components (i.e., oxygen and plankton). Arguably, in the marine environment, the most common and the most powerful of these processes is the turbulent mixing. Having taken its effect in the form of the turbulent diffusion [15], System (4)-(6) turns into the following: ∂u(x, y, t) ∂v(x, y, t) (in dimensionless units) where the turbulent diffusion coefficient is equal to one due to our choice of dimensionless spatial coordinates [12]. Plankton system is affected by several environmental factors. These factors can be physical, i.e., light, temperature, salinity, turbulence etc., or biological such as predation and availability of nutrients [16][17][18][19]. Keeping in mind the possibility to link the change in the dynamics to the effect of the global warming (that can affect the rate of oxygen production) [10,12], we consider the properties of the emerging spatiotemporal patterns subject to the value of the net oxygen production rate A, which we therefore regard as a controlling parameter. In what follows, pattern formation is studied by means of extensive numerical simulations for a different choice of the initial conditions and for a sequence of values of A approaching (and eventually over-passing) the critical value when the sustainable dynamics becomes impossible.

Simulation Results: Pattern Formation
Equations (7)-(9) were solved numerically by the finite-differences method in a rectangular domain with 0 < x < L x and 0 < y < L y . At the domain boundary, the no-flux Neumann boundary conditions were imposed. Simple explicit scheme was used. It was checked that the choice of the spatial and temporal mesh steps is sufficiently small to avoid any essential numerical artifacts. For simulations, we chose parameter values as h 1 = 0.7, h 2 = 1, h 3 = 1, h 4 = 1, η = 0.01, B = 1.8, γ = 1.2, β = 0.7, µ = 0.1, h = 0.1 and σ = 0.1, and vary A over a certain range. Our choice of parameters is based on the previous studies of the corresponding plankton-oxygen model with one spatial dimension [10,12] where it was shown that the above parameter set ensures (up to some restrictions on possible values of A) the sustainable functioning of the system.
Simulations are performed for a few different types of initial conditions. Our first choice of the initial conditions describes spatially homogeneous distributions for both oxygen and phytoplankton at their steady states whilst the zooplankton density is a linear function of space with a very small gradient: where the values used in simulations are 2 = 3 × 10 −5 and 3 = 6 × 10 −5 . We mention here that, based on the properties of some other well-known and similar models (e.g., see [20,21]), the details of the initial species distribution can affect the system's dynamics at small time but they are unlikely to affect the dynamics at large time, after the transients associated with the initial conditions disappear. Some snapshots of the system dynamics for the initial conditions given by Equations (10)-(12) are shown in Figures 2 and 3. The distributions of oxygen, phyto-and zooplankton exhibit qualitatively similar behavior, hence only the oxygen concentration is shown. Figure 2 (obtained for A = 2.09) shows that the evolution of the initial conditions first leads to the formation of stripes. Eventually the stripes are destroyed by the perturbation induced by the boundary resulting, in the large-time limit, in an irregular patchy structure.
Computation 2018, 6, 59 5 of 14 the initial species distribution can affect the system's dynamics at small time but they are unlikely to affect the dynamics at large time, after the transients associated with the initial conditions disappear. Some snapshots of the system dynamics for the initial conditions given by Equations (10)-(12) are shown in Figures 2 and 3. The distributions of oxygen, phyto-and zooplankton exhibit qualitatively similar behavior, hence only the oxygen concentration is shown. Figure 2 (obtained for A = 2.09) shows that the evolution of the initial conditions first leads to the formation of stripes. Eventually the stripes are destroyed by the perturbation induced by the boundary resulting, in the large-time limit, in an irregular patchy structure. The dynamics become different for a larger value of A. Figure 3 shows snapshots of the oxygen concentration arising from the same initial conditions Equations (10)-(12) but for A = 2.3. In this case, the stripes that emerge at the early stage of dynamics are not being destroyed but evolve to spirals.
The spiral pattern appears to be sustainable (which was checked in simulations running up to t = 6000) and is not evolving to the patchy structure. We mention here that for this value of A the nonspatial system is unsustainable regardless of the initial conditions (resulting in the oxygen depletion and plankton extinction) [12].
We now consider a different initial condition where only the oxygen concentration is uniform but both the phyto-and zooplankton densities are described by slowly changing functions: where 1 = 2 × 10 −7 , 2 = 6 × 10 −7 and 3 = 3 × 10 −5 . For these initial conditions, the formation of irregular patchy structure in the large-time limit is shown in Figure 4 (obtained for A = 2.06). Since in this case the initial conditions are lacking planar symmetry, no stripes are seen at the intermediate time. The dynamics become different for a larger value of A. Figure 3 shows snapshots of the oxygen concentration arising from the same initial conditions Equations (10)-(12) but for A = 2.3. In this case, the stripes that emerge at the early stage of dynamics are not being destroyed but evolve to spirals.
The spiral pattern appears to be sustainable (which was checked in simulations running up to t = 6000) and is not evolving to the patchy structure. We mention here that for this value of A the nonspatial system is unsustainable regardless of the initial conditions (resulting in the oxygen depletion and plankton extinction) [12].
We now consider a different initial condition where only the oxygen concentration is uniform but both the phyto-and zooplankton densities are described by slowly changing functions: v where 1 = 2 × 10 −7 , 2 = 6 × 10 −7 and 3 = 3 × 10 −5 . For these initial conditions, the formation of irregular patchy structure in the large-time limit is shown in Figure 4 (obtained for A = 2.06). Since in this case the initial conditions are lacking planar symmetry, no stripes are seen at the intermediate time.  However, for a larger value of A the patterns become different showing the same tendency as above, i.e., they become more regular and exhibit a clear spiral structure. A typical example is shown in Figure 5 obtained for A = 2.4. Similar dynamics are observed for some larger values of A, e.g., see Figure 6 obtained for A = 2.7 and the same initial conditions given by Equations (13)- (15).  However, for a larger value of A the patterns become different showing the same tendency as above, i.e., they become more regular and exhibit a clear spiral structure. A typical example is shown in Figure 5 obtained for A = 2.4. Similar dynamics are observed for some larger values of A, e.g., see Figure 6 obtained for A = 2.7 and the same initial conditions given by Equations (13)- (15). However, for a larger value of A the patterns become different showing the same tendency as above, i.e., they become more regular and exhibit a clear spiral structure. A typical example is shown in Figure 5 obtained for A = 2.4. Similar dynamics are observed for some larger values of A, e.g., see Figure 6 obtained for A = 2.7 and the same initial conditions given by Equations (13)- (15).
The spiral pattern appears to be stable and sustainable as was checked in simulations running up to t = 6000. Formation of the spirals at large time is preceded by the transient snake-like pattern (cf. the top-right panel in Figures 5 and 6). We want to emphasize that, for values A ≥ 2.4, the corresponding 1D system is not sustainable and the evolution of initial distributions promptly leads to the oxygen Computation 2018, 6, 59 7 of 14 depletion and the plankton extinction [12]. (Here we also recall that the corresponding nonspatial system becomes unsustainable already for A ≥ 2.055, see [12].) Therefore, the higher persistence of the system observed in the 2D case should be attributed solely to the effect of the pattern formation. The spiral pattern appears to be stable and sustainable as was checked in simulations running up to t = 6000. Formation of the spirals at large time is preceded by the transient snake-like pattern (cf. the top-right panel in Figures 5 and 6). We want to emphasize that, for values A ≥ 2.4, the corresponding 1D system is not sustainable and the evolution of initial distributions promptly leads to the oxygen depletion and the plankton extinction [12]. (Here we also recall that the corresponding nonspatial system becomes unsustainable already for A ≥ 2.055, see [12].) Therefore, the higher persistence of the system observed in the 2D case should be attributed solely to the effect of the pattern formation.  depletion and the plankton extinction [12]. (Here we also recall that the corresponding nonspatial system becomes unsustainable already for A ≥ 2.055, see [12].) Therefore, the higher persistence of the system observed in the 2D case should be attributed solely to the effect of the pattern formation.  Now we are going to check whether the results obtained above are sensitive to the size of the domain and to a further modification of the initial distributions. We therefore perform simulations in a larger domain 900 × 900 for the following initial conditions: where 1 , 2 and 3 have the same values as in (13)- (15). Figure 7 shows snapshots of the oxygen concentration obtained for A = 2.04 and the initial conditions given by Equations (16)- (18). We observe that, in the large time run, the system eventually evolves to an irregular patchy distribution similar to the ones shown in the bottom-right panel in Figures 2 and 4. There are, however, some differences: in the case shown in Figure 7 (middle row), the irregular patchy pattern is preceded by the formation of a double spiral or a "mushroom-like" structure. For a somewhat larger value of A, the dynamics exhibits similar features (see Figure 8 obtained for A = 2.07) except for the relatively minor difference that the spirals are of a different shape and that for the intermediate time the spatial distribution is becoming more prominently heterogeneous, cf. the dark-blue areas in the middle and bottom rows of Figure 8 and the range of values on the color bar.   An increase in A results in a change in the emerging pattern similar to what was observed in the smaller domain (cf. Figures 2 and 3), i.e., the evolution of the initial conditions leads to the formation of a clear spiral pattern, see Figure 9 obtained for A = 2.4 and the initial conditions given by Equations (13)- (15). The pattern is stable and sustainable and does not evolve into an irregular or patchy distribution; this was checked in simulations running up to t = 15,000. An increase in A results in a change in the emerging pattern similar to what was observed in the smaller domain (cf. Figures 2 and 3), i.e., the evolution of the initial conditions leads to the formation of a clear spiral pattern, see Figure 9 obtained for A = 2.4 and the initial conditions given by Equations (13)- (15). The pattern is stable and sustainable and does not evolve into an irregular or patchy distribution; this was checked in simulations running up to t = 15,000.
A further increase in A to somewhat larger values (up to about A ≈ 3.5) does not change the dynamics in a significant way but can make the spirals or mushroom-like patterns be seen more clearly at already smaller time. Figure 10 shows snapshots of the oxygen concentration obtained for A = 2.5. The evolution of the initial conditions first leads to a zigzag, mushroom-like structure (already at t ∼ 100) from which a clear double spiral eventually develops. The double spiral is stable and sustainable and is not further evolving to a distribution with different properties (e.g., patchy) as was checked in long-run simulations.
Once again, we notice that the persistent spatial patterns shown in Figures 9 and 10 exist in the parameter range where the system becomes unsustainable in the corresponding 1D case (i.e., for A ≥ 2.4). We therefore consistently observe that the dynamics of the 2D system is sustainable in a broader parameter range. However, a much larger increase in A (approximately up to 3.5 or larger) will make the 2D system non-viable too, resulting in plankton extinction and oxygen depletion. A further increase in A to somewhat larger values (up to about A ≈ 3.5) does not change the dynamics in a significant way but can make the spirals or mushroom-like patterns be seen more clearly at already smaller time. Figure 10 shows snapshots of the oxygen concentration obtained for A = 2.5. The evolution of the initial conditions first leads to a zigzag, mushroom-like structure (already at t ∼ 100) from which a clear double spiral eventually develops. The double spiral is stable and sustainable and is not further evolving to a distribution with different properties (e.g., patchy) as was checked in long-run simulations. A further increase in A to somewhat larger values (up to about A ≈ 3.5) does not change the dynamics in a significant way but can make the spirals or mushroom-like patterns be seen more clearly at already smaller time. Figure 10 shows snapshots of the oxygen concentration obtained for A = 2.5. The evolution of the initial conditions first leads to a zigzag, mushroom-like structure (already at t ∼ 100) from which a clear double spiral eventually develops. The double spiral is stable and sustainable and is not further evolving to a distribution with different properties (e.g., patchy) as was checked in long-run simulations.  We also notice that in the parameter range where the dynamics is only sustainable in a higher-dimensional system (e.g., 2D vs. 1D), the dynamics becomes more sensitive to the details of the initial conditions. For a given parameter set, whilst one choice of the initial distributions would evolve to a sustainable spiral pattern, another choice might eventually lead to extinction. One example We also notice that in the parameter range where the dynamics is only sustainable in a higher-dimensional system (e.g., 2D vs. 1D), the dynamics becomes more sensitive to the details of the initial conditions. For a given parameter set, whilst one choice of the initial distributions would evolve to a sustainable spiral pattern, another choice might eventually lead to extinction. One example is shown in Figure 11. The evolution of the initial distributions leads to the formation of a ring-like elliptic shape. The shape then grows in size hence 'propagating' through the space towards the domain boundaries. Outside and inside the ring, the oxygen concentration and the plankton densities are close to zero. Eventually the ring collides with the boundary and disappears, thus resulting in the system collapse (extinction and depletion).
Computation 2018, 6, 59 11 of 14 is shown in Figure 11. The evolution of the initial distributions leads to the formation of a ring-like elliptic shape. The shape then grows in size hence 'propagating' through the space towards the domain boundaries. Outside and inside the ring, the oxygen concentration and the plankton densities are close to zero. Eventually the ring collides with the boundary and disappears, thus resulting in the system collapse (extinction and depletion).

Discussion and Concluding Remarks
Dynamics of the dissolved oxygen in the ocean has been attracting considerable attention over the last decade [8,9,[22][23][24]. The existing observations point out at a gradual decrease of the average oxygen level across the globe, with the areas of extreme conditions (a severe lack of oxygen resulting in anoxia and mass mortality of the marine fauna) being fast growing in size. Although there is a consensus that the oxygen decrease is a result of global climate change, in particular global warming, the specific reasons and mechanisms of the decrease remain controversial. The decrease was previously linked to purely physical mechanisms such as the decrease of oxygen solvability in warmer waters [23,24]. More recently, however, it was shown that global warming can significantly affect the amount of dissolved oxygen also through a biological mechanism as it can disrupt the oxygen production by phytoplankton potentially leading to global ocean anoxia [10,12,25]. Moreover, since it is estimated that about 70% of the atmospheric oxygen is produced in the ocean, such a global disruption of oxygen production would eventually lead to a significant decrease of the oxygen concentration in the air and could hence lead to mass mortality of animals and humans.
In our previous work, the above issue was addressed theoretically using a novel mathematical model of the coupled plankton-oxygen dynamics and extensive computer simulations [11,12]. However, the specific focus was on the dynamics of non-spatial ("well-mixed") system and the system with one spatial dimension. For the nonspatial system, it was shown that the sustainable dynamics is Note that all parameters are the same as in Figure 5, but the initial conditions are slightly different, i.e., c(x, y, 0) = 0.4723, u(x, y, 0) = 0.3612 − 2 × 10 −7 (x − 180)(x − 420) − 6 × 10 −7 (y − 90)(y − 210) and v(x, y, 0) = 0.1554 − 6 × 10 −7 (x − 250) − 3 × 10 −5 (y − 250). Contrary to the case shown in Figure 5, now the evolution of the initial conditions results in oxygen depletion and species extinction.

Discussion and Concluding Remarks
Dynamics of the dissolved oxygen in the ocean has been attracting considerable attention over the last decade [8,9,[22][23][24]. The existing observations point out at a gradual decrease of the average oxygen level across the globe, with the areas of extreme conditions (a severe lack of oxygen resulting in anoxia and mass mortality of the marine fauna) being fast growing in size. Although there is a consensus that the oxygen decrease is a result of global climate change, in particular global warming, the specific reasons and mechanisms of the decrease remain controversial. The decrease was previously linked to purely physical mechanisms such as the decrease of oxygen solvability in warmer waters [23,24]. More recently, however, it was shown that global warming can significantly affect the amount of dissolved oxygen also through a biological mechanism as it can disrupt the oxygen production by phytoplankton potentially leading to global ocean anoxia [10,12,25]. Moreover, since it is estimated that about 70% of the atmospheric oxygen is produced in the ocean, such a global disruption of oxygen production would eventually lead to a significant decrease of the oxygen concentration in the air and could hence lead to mass mortality of animals and humans.
In our previous work, the above issue was addressed theoretically using a novel mathematical model of the coupled plankton-oxygen dynamics and extensive computer simulations [11,12]. However, the specific focus was on the dynamics of non-spatial ("well-mixed") system and the system with one spatial dimension. For the nonspatial system, it was shown that the sustainable dynamics is only possible when the oxygen production rate (i.e., A) is in an intermediate range, i.e., not too high and not too low [10,12]. Interestingly, it was shown that the corresponding 1D spatial system exhibits similar properties but the parameter range of sustainable dynamics is somewhat broader. Therefore, the question arises as to whether the sustainability range may be even broader in a more realistic 2D system (i.e., system with two spatial dimensions) or perhaps the oxygen depletion catastrophe would not happen at all. Indeed, there are many examples in theoretical ecology where the 2D system is sustainable whilst the corresponding 1D or nonspatial systems are not; e.g., see Chapter 12 in [20].
In this paper, we have considered the above question using the 2D spatially-explicit model of the plankton-oxygen dynamics; see Equations (7)- (9). Our main focus is on the pattern formation as it is well known from other population dynamics models that pattern formation can enhance system's sustainability [20]. Since our ultimate goal is to put our results into the context of global warming (cf. [10,12]), we use A as the controlling parameter as there is significant biological evidence that the rate of net oxygen production depends on the water temperature [26][27][28][29]. Using extensive numerical simulations, we have shown that the model exhibits rich spatiotemporal dynamics. In the parameter range where the nonspatial system either possesses a stable limit cycle or becomes non-sustainable when the cycle is destroyed in a non-local bifurcation (see [10] for details), the evolution of the initial conditions leads to the formation of spatiotemporal patterns. In particular, we have observed (for somewhat different initial conditions) irregular spatiotemporal patterns, simple spirals, double spirals or "mushroom-like" structures, and "snake-like" structures. Here we recall that the spatial distribution of the dissolved oxygen typically shows a patchy structure [8,9], although the quality of the available data (e.g., the spatial resolution) does not make it possible to look into the details of that structure. Whilst quantitative comparison between theory and data lies beyond the scope of this paper, we conclude that our results are in qualitative agreement with at least some of the field data. We also mention here that spiral-like and mushroom-like spatial plankton patterns are frequently observed in the ocean [30,31]. Mushroom structures are usually linked to particular hydrodynamical conditions (e.g., the meeting of two masses of water, one cold and the other warm, cf. [30]); the results of our study apparently suggest an alternative explanation where the formation of such structures mostly results from biological interactions.
Interestingly, the type of the pattern depends on A. Whilst a value of A taken in the middle of the sustainable range results in a nearly-uniform oxygen and plankton distribution, for larger values of A normally lead to spatially irregular, strongly heterogeneous, large-amplitude patterns. Along with a further increase in A, patterns tend to become more regular. For a sufficiently large value of A, the depletion of oxygen and the plankton extinction (which can be preceded by long living spiral patterns) is the only issue, similar to what was observed in the corresponding 1D and nonspatial models. The value of A for which the extinction is observed appears to be slightly larger than it is in the 1D case.
We also notice that for larger values of A-i.e., the values that are close to the sustainable parameter range-the spiral or double-spiral patterns appear to be typical. This seems to suggest that the emergence of spirals in the spatial distributions of the oxygen concentration can be regarded as an early warning signal of the approaching tipping point.
Our study leaves a few open questions. In this paper we considered the two spatial dimensions with no spatial structure, which arguably corresponds to the horizontal coordinates. This seems to be a natural extension of the 1D case. An interesting alternative could be a 2D model describing the vertical distribution of oxygen and plankton, i.e., where one spatial coordinate would correspond to the depth. Such a model should include an explicit dependence on the vertical coordinate in order to account for the ocean stratification. One can expect that such a system could have rather different properties.
Note that in this paper our goal is to reveal some generic scenarios of pattern formation as it helps to understand the limits of the system's sustainability and also to identify possible early warning signals of the approaching tipping point. For this purpose, dimensionless variables and parameters are sufficient. However, in order to make the study more quantitative and to eventually compare our predictions to field or laboratory data and, ultimately, to increase the model's predictive power, the results should be presented using the original, dimensional physical/biological units.
Technically, it is easy to do using the definition of the dimensionless variables-as soon as all the original parameters are known. The latter, however, is a challenging task and requires much more than routine work as the value of most of the parameters is known with a low precision, sometimes only up to an order of magnitude. Making a careful and sensible analysis of the available empirical data to translate our results to the real, dimensional units lies beyond the scope of this paper and but will become a focus of future work.
Another open issue is the description of the turbulent transport. Although the turbulent diffusion framework is a widely used approach, it has its limitations, especially in a multiscale, fully developed turbulent flow. One alternative framework is based on the theory of chaotic flows [32]. Several studies on the population dynamics in marine ecosystems [33,34] used this framework and showed that its predictions can be different from those obtained using the turbulent diffusion approach. Extending our model of the coupled plankton-oxygen dynamics to better describe the properties of the turbulent mixing is a challenge for future work.
Author Contributions: S.P. designed the research. Y.S. wrote the numerical code and performed simulations. S.P. and Y.S. wrote the manuscript.