The Inﬂuence of Reduced Light on Competing Zooplanktonic Populations – Inference from Mathematical Modeling and Field Data

A diffusion model for the competition of three zooplankton populations is developed, that is based on earlier ideas on reduced phytoplankton abundance at large depths in the ocean, due to light absorption in the higher water layers and shading by other upper residing microorganisms. Equilibria of the diffusion-free system are analysed. A sensitivity analysis in terms of some population interaction parameters that are difﬁcult to assess is performed. Simulations of the space-dependent system reveal the changes in the populations abundances in the various water layers, showing their expected reduction with increasing depths.


Introduction
The knowledge of the planktonic organisms to understand the dynamics of the marine ecosystem is highly important, as well as the study of the zooplanktonic crustacean copepod community is crucial for monitoring changes in marine ecosystems, [5,24]. Indeed within the pelagic realm, planktonic copepods are among the most abundant and widely distributed mesozooplanktonic communities and main primary consumers, making up about one third of the mean plankton biomass and half of the mean zooplankton abundance in the marine environment, [16] Moreover marine copepods have an essential role in the function of marine trophic webs and they are central to the productivity and biogeochemistry of marine ecosystem. Furthermore copepods respond sensitively to environmental changes and have been recognized as indicators of anthropogenic stressors, [32,14,13,3,27]. The perception of the trophic role of copepods has changed during recent decades. In the past, they were viewed as purely herbivorous, which acted as a link between primary production and planktivorous fish. Copepod food selection does not only depend on food particle size, [28]. Indeed it has been shown that copepods can strongly select between similarly sized food particles because of chemical properties, [9], an ability which is lacking in true zooplanktonic filter feeders such as tunicates and cladocerans. Copepods are thought to use different feeding modes depending on their species-specific mouthparts, particularly in the morphology of the mandibular gnathobases and the movements that are possible with them [26,22,21,4]. More generally, is widely known that zooplankton mainly feeds on phytoplankton, hence the phytoplankton plays an important role in the study of zooplankton dynamics.

559
Copepods can potentially obtain food from any known stock of organic matter in dissolved or particulate form; their feeding appendages and behaviour allow them to capture particles, phytoplankton or detritus, of just a few microns in size or to attack living zooplankton organisms, such as chaetognaths, medusae or other copepods, [23,6]. According to their main dietary is possible to recognize three distinct trophic groups in mesozooplanktonic copepods: herbivores, which feed mainly unicellular algae and diatoms with silicified frustules, omnivores able to feed both phytoplankton and micro-mesozooplankton and carnivorous able to prey herbivorous and omnivorous copepods. In this study we assume that the herbivores and omnivores are prey only for the carnivorous copepods. Omnivores are not able to prey herbivorous copepods due their sizes. Among the environmental parameters, the availability and intensity of sunlight plays a pivotal role. Indeed, the variation of light dependent to the season, deep and water turbidity may have indirect effects on the feeding strategies, food availability, competition and efficiency in copepods, with ecological implications for the functioning of the whole ecosystem, [31,30,15]. Note that this seasonality effect is incorporated in our model through the light attenuation function. Indeed, the latter has been obtained from data sampled in each different season, as it will be better clarified when explaining the model. For zooplankton grazers, including copepods, light offers a behavioural dilemma that is potentially present at both small and large scales. Copepods rely on phytoplankton that in turn rely on light for photosynthesis. In the same time copepods need to avoid being detected by visually feeding fish that see better in light environments. That is especially the case of the Diel Vertical Migration, whereby zooplankton descend to depth during daylight hours and ascend again to feed during darkness. The intensity of migrations is mainly related to light intensity and predation risk, [12,18]. Direct support for this notion comes from experiments in which the amplitude of migration varied in proportion to predator abundance and predation pressure, [7]. Furthermore, Transparency Regulator Hypothesis (TRH), developed in [33] combines biotic (food resources and visual predation) and abiotic (temperature and ultraviolet radiation) factors regulating presence and abundance of zooplankton through the water column.

Study area and sampling method
The sampling site is located in the neritic waters off the Tuscan coast in the Southern Ligurian Sea at the North East border of the Western Mediterranean and is connected to the southern basin (Tyrrhenian Sea) across the Corsica Channel. The major large-scale feature of the water dynamic of the Ligurian Sea is a cyclonic circulation that is active all year round, but more intense in winter than in summer, involving both deep and surface layers [1]. Climatic forcing can greatly change the intensity of fluxes, but the general pattern can be considered permanent, [17,19]. The sampling station (4328'10"N, 1001'55"E) is characterized by the large extension of the continental shelf and limited depth (100 m), even at remarkable distances from the coast (12.5 nautical miles), [4]. The sampling area represents the water column for the models. The entire water column is sampled with a WP-2 standard net (mesh size 200 m and diameter 57 cm) from the surface to 100 m of depth. The samples were fixed in 4% neutralized formaldehyde in sea water buffered with borax and kept in the dark, [8]. Due to the complexity of phytoplankton dynamics analisys, in addition to those of zooplankton, in this work phytoplankton is taken as an unlimited resource.

The model without diffusion
We want to model the copepods interactions via the mathematical model that is introduced below. The model is built along the ideas exposed in [2] which date back also to [20,11]. In practice, a term denoting the light absorption and the decrease of its availability in the layers of the ocean is introduced. In [2] this is explicitly accounted for the phytoplankton reproduction. Observe that the light attenuation is modeled by a least squares fit on four sets of real sampled data, that have been taken in the four seasons. Thus the parameter ν indicates the decrease in luminosity at the various depths, but it has to be understood to have a different value in each season. Here, we do not explicitly model phytoplankton, but we take anyway into consideration the fact that it becomes less abundant in the deeper the layer in the ocean, with consequences for zooplankton feeding. Let H, C and O respectively denote the herbivores, carnivores and omnivores populations at the depth x, which assumed to be fixed in the non-spatial model. Thus this model describes the ecosystem dynamics at every and each different depth in the water column.
The copepods dynamics is given by the following system in which all the parameters are nonnegative: In this system the resource that is at the bottom of the trophic chain, phytoplankton, is present and assumed to be unlimited, but it is not explicitly modeled. The first term denotes, as commented above, the light attenuation with the ocean depth and it takes different values in each season, these being obtained by a linear regression on actual sampled field data, that are used in the simulations carried out later. This models the fact that light absorption is a limiting factor for photosynthesis, a fact that becomes important for the phytoplankton sustainable presence. The herbivore copepods growth rates therefore completely depend on it, while the omnivores only partially, but are anyway substantially affected as well.
The first term in the first two equations therefore contains the exponential factor in which x denotes the water depth and ν the attenuation coefficient. In the first equation the next two terms denote competition among copepods and herbivores, herbivores natural mortality. The final Holling type II term represents the predation of carnivores over herbivores.
The second equation contains the dynamics of the omnivores, that can feed both on phytoplankton as mentioned above, first term, and have alternative resources on which they feed at rate p, represented by the second term, their natural mortality rate being q. Carnivores in addition feed on omnivores, this being represented by the Holling type II term. The last two terms indicate losses for the class of omnivores, having a negative sign. They are not due to direct predation, but rather to indirect competition for resources among omnivores and the other two zooplankton populations, namely with carnivores at rate j and herbivores at rate c. This is competition, as these terms are matched by other negative terms in the equations for H and C, specifically those with coefficients a and g. In view of the fact they the populations abundant, we model these actions by simple mass action terms.
The third equation shows the carnivores dynamics. The first and the last term represent the logistic carnivores growth due to the presence of types of zooplankton other than copepods, at rate k, and intraspecific competition at rate n. The next two terms represent the gain obtained by predation on both herbivores and omnivores, scaled via respective conversion factors e and v. The fourth term denotes the carnivores competition with the omnivores at rate g.

Analysis
The system can settle to several different states: but note at first that the one in which three copepods populations are absent and the ecosystem disappears cannot be attained as it is unstable. Similarly the carnivore-free point is not stably attainable.
The steady states that can be achieved under different parametric conditions are as follows: the carnivores-only point E 3 = (0, 0, C 3 ), the omnivores-free point E 5 = (H 5 , 0, C 5 ), the herbivores-free point E 6 = (0, O 6 , C 6 ) and coexistence E 7 = (H 7 , O 7 , C 7 ). The feasibility of these equilibria is studied in Appendix A, while their stability results are reported in Appendix B. At E 3 the population level is C 3 = kn −1 ; it is stable for The omnivores-free situation E 5 can give rise to one or two points. The feasibility condition for exactly one equilibrium is whose population values H 5+ and C 5+ are given in (6) below. This equilibrium does not occur at all if instead Two such equilibria could also exist if In such case the population levels are Stability is achieved under the conditions In the case of the herbivores-free point E 6 , A unique equilibrium is guaranteed for while two such equilibria exist if For stability, the following conditions must be satisfied L < aO 6 + bC 6 , n > wz Equilibrium points Feasibility conditions Stability conditions .
,  Finally, ecosystem coexistence E 7 will be investigated numerically. We report now our findings on ecosystem bifurcations. It is observed that near each and every one of these equilibria, the system undergoes possible local bifurcations. These results are formally proved in Theorems 1-6 in the Appendix. In summary, we find the following situations.
At E 3 , taking r as bifurcation parameter, when it crosses the critical value there is a transcritical bifurcation, but neither saddle-node nor pitchfork bifurcations. Taking instead the bifurcation parameter k, the system exhibits a pitchfork bifurcation but neither transcritical nor saddle-node bifurcations, for This is a semicritical situation, since the parameter k by definition must of course be nonnegative. At E 5 we find again in terms of the bifurcating parameter r that there no local bifurcations, even though the parameter attains the critical value At E 6 again in terms of r there is a transcritical bifurcation but neither saddle-node nor pitchfork bifurcations for the critical value In terms of the bifurcation parameter λ = k for the critical value a saddle-node bifurcation arises, but neither transcritical bifurcation nor pitchfork bifurcations occur.
Finally for E 7 there is also a saddle-node bifurcation and neither transcritical bifurcation nor pitchfork bifurcations when µ = r attains the value

Simulations
We consider one dimensional the diffusion model related with (1), x being the vertical depth, measured from the upper ocean surface downwards, where We have discretized it with a suitable numerical forward finite difference scheme, taking care to choose the time advancement step small enough so as to satisfy the von Neumann stability condition. We chose a reference set of parameter values, the same for all the simulations. In view of the fact that the spatial system does not settle to an equilibrium, we ran the simulations for 10 days, for all the situations. At that time, we then saved and plotted the zooplankton population distributions along the vertical column, at the same depths, arbitrarily chosen at 1m, 10m, 15m, 25m, 50m and 75m. In addition we calculated and plotted the average value of each population over the whole water column. We repeated this procedure by suitably secting pairs of parameters that are difficult or impossible to assess in the laboratory or in the field, mostly the interspecific competition coefficients, and we let them vary in a reasonable range, to establish the sensitivity of the solution at the specified time with respect to these parameters. The plots contain surfaces as function of two parameters, as stated above. The height of the surface at each point in the domain denotes the population level of the system for the combination of the parameter values that are the coordinates of that point. In this way it is possible to see if and how the population level changes with respect to certain perturbations in the parameter values.
This procedure has been repeatead in each case also for the four seasons. We now illustrate a sample of the results.

Specific results for the parameter space a − z
We consider now the case of changes in the parameters a ∈ [0.01, 0.1], namely the damage felt by the herbivores in the competition of the omnivores and z ∈ [0.1, 0.7], the pressure exerted on omnivores by the carnivores. Figure  1 shows the total population values in each season.
Focusing first on the winter season, concentrating on the whole populations over the entire water column, we observe that the herbivores are by far the most abundant population, at levels reaching 4.5 × 10 4 , followed by omnivores, at about one order of magnitude lower, 1000 and then by carnivores at 120. The shape of the distribution shows that the herbivores peak is located for low values of a, independently of the values of z. Instead, the omnivores have a peak for higher values of a and lower values of z, while carnivores still exhibit a different shape, with high values concentrating over larger z and almost independent of a, but with a maximum at an intermediate value of a.
With respect to winter, in the spring all the populations increase and their surfaces retain the same qualitative behavior. The herbivore population attains a value higher, 5 × 10 5 the omnivores show a maximum at 2500; the carnivores are found at a maximal value 140.
In the summer, with respect to their values in the spring, again the shapes of the surfaces are very much alike. The maximal values become even larger, for the herbivores by one order of magnitude, 1.2 6 ; the omnivores attain the level 3500, the carnivores are roughly at the same level of the spring.
In the fall a reverse situation with respect to what happens in the summer is observed. Indeed, comparing values with the summer season, the herbivores drop to 7 × 10 5 , the omnivores also drop to 3500, while the carnivores experience a slight increase in some parts of the parameter range, but with a maximum that does not change, at 140.
The transition from fall to winter, shows another reduction of all the populations, herbivores are halved, omnivores drop by a third, carnivores also are reduced by roughly 20%. Now we examine each depth, descending along the water column.

Winter
In the winter at depth 1m, both the population values and the shapes as functions of the two parameters follow the trends of each whole population along the water column, this feature being repeatead at every depth, except possibly the very deep layer, in which however the the values are low anyway and such changes do not matter much. Herbivores dominate at value 4000, followed by omnivores at 60 and lastly by carnivores with a peak at 2.5. At 10m the herbivores and omnivores fall by about a third to respective maxima of 2000 and 40, carnivores are essentially unchanged. This trend continues at 15m, herbivores are found at 1200, omnivores at 30, carnivores unchanged. At 25m, herbivores are found at 250, omnivores drop to 12, carnivores again do not change. At 50m, herbivores drop by two orders of magnitude to 3.5, omnivores are also found in the single digit range at 2.5, carnivores are halved, attaining the value of 1.4. At 75 herbivores are reduced by a further order of magnitude the lowest rank among the populations at 0.25, omnivores fall to 0.45, carnivores thrive around 0.35.

Spring
In the spring at depth 1m the shapes of the three surfaces for the populations are again similar to the respective ones of the total populations, herbivores and carnivores having higher values for low a, independently of z, while omnivores do the reverse. They attain maxima respectively at 3.5 × 10 4 , 140, 2.5. These values drop at 10m for the herbivores and omnivores to 2.5 × 10 4 and 120 with carnivores unchanged. At 15m the trend follows with herbivores and omnivores at 1.4 × 10 4 and 80, carnivores still unaffected. At 25m the herbivores rebound to values 3000 while omnivores continue to decrease at 35 and carnivores are not affected. A marked decrease is found at 50m for the whole ecosystem, with populations respectively attaining 30, 4 and 1.4 and continues at depth 75m with values now settling to 0.4, 0.5 and 0.35.

Summer
In the summer again in general the surface shapes follow the whole population distributions. Figure  2 shows the abundance of each zooplankton population at the selected depths. At 1m herbivores have a peak at 8 × 10 4 , omnivores at 180, carnivores at 2.5. At 10m, we find the herbivores slightly reduced at 6

Fall
In the fall the shapes of the surfaces at each depth follow the one of the average distribution. At 1m herbivores are the most abundant population, but much reduced with respect to the sumer. Theis peak is found at 5000, followed by omnivores at 70 and carnivores at 2.5. The first two populations decrease at 10m with respective peaks at 3500 and 50, while carnivores remain at 2.5. This trend continues at 15m with 2000 and 40, at 25m with 500 and 40, with carnivores always unchanged. At 50m we respectively find 9, 3.5 and 1.4, and finally at 75m the values are 0.4, 0.5, 0.35.

Specific results for the parameter space j − z
We consider now the case of the parameters for the competition j of carnivores over omnivores and the direct predation of carnivores on omnivores z. For the whole populations, see Figure 3, in the winter the omnivores are the predominant population, with maximal values at 700, herbivores are situated at 140 and carnivores at 120. The parameter j seems to have less influence for the latter and for omnivores. Carnivores are depressed for low values of z, while omnivores behave in the opposite way, in particular they tend to increase or decrease when both parameters are simultaneously reduced or become larger. In the opposite way, herbivores increase when both parameters have higher values and drop when anyone of the two decreases. In the spring there is an increase in all the populations, herbivores peak to 400, omnivores at 1200, carnivores at 140. For the summer the peaks are found even at higher values as expected, respectively at 900, 1400, carnivores again at 140. Finally in the fall there is a drop, except for carnivores, with respective values at 120 and 900, followed again with another decrease in the winter for omnivores and carnivores, by comparison with the values mentioned above, while herbivores instead experience a slight increase, perhaps due to less competition by the other populations. For the various depths, independently of the season, the population surfaces behave like their average counterparts described above. We now examine for each season the differences observed in each depth.

Winter
In the winter, at 1m the omnivores dominate with a maximum at 35, herbivores are 10, decreasing respectively to 24 and 6, 18 and 3.5 and 1 and 8 at 10m, 15m and 25m, while carnivores are always found at level 2.5. In the deeper layers, 50m and 75m, all populations drop further, omnivores at 2 and 0.35, herbivores at 0.22 and 0.08, carnivores 1.4 and 0.3.

Spring
The increase of the ecosystem populations, with respect to the winter season, at 1m leads to the values 60 for the omnivores, 25 for herbivores and 2.5 for carnivores. At 10m the first two population values respectively reduce to 50, 18 and carnivores are unaffected. At 15m the trend continues with values 35 and 12, followed by 18 and 4 at 25m, with no changes in the carnivores. At 50m and 75m we find omnivores at 3.5 and 0.5, herbivores at 0.5 and 0.13, carnivores at 1.4 and 0.35.

Summer
In the summer, at 1m the omnivores are at 70, herbivores at 60, decreasing respectively to 50 and 45, 40 and 25 and 20 and 8 respectively at the depths 10m, 15m and 25m, while carnivores remain always at 2.5. At 50m omnivores are at 3.5, herbivores at 0.55, carnivores at 1.4, and further down at 75m, these values respectively become 0.5 and 0.13 and 0.35.

Fall
In the fall, the same configuration occurs, although with respect to the summer the populations decrease. At 1m the omnivores are at 40, herbivores 8 and these values decrease respectively to 30

Specific results for the parameter space b − g
We consider now the case of changes in the parameters b, namely the damage felt by the herbivores in the competition of the carnivores and g, the pressure exerted on carnivores by the omnivores. Concentrating on the whole populations over the entire water column, see Figure 3, we focus at first on the winter season. We observe that the omnivores are the largest population by one order of magnitude, with the peak at 550, the herbivores and carnivores one order of magnitude lower, respectively with maximum values about 75 and 40.
In the spring there is the usual increase in the numbers, but the omnivores remain the most abundant population, at 1070, herbivores reach 95, carnivores 45. A further increase is observed in the summer with respective values attaining 1130, 130 and carnivores remaining at 45. The expected decline follows in the fall, with the maximal populations values being 745 for omnivores, 75 for herbivores and 40 for carnivores, the least affected ones.
The surface shapes indicate three different trends. Namely, for low values of b and g, herbivores peak, markedly decreasing especially for larger values of b; omnivores increase with larger g and are not much affected by changes in b; carnivores have maximal values for high b and low g. Now for every season we examine each depth in detail.

Winter
In the winter at depth 1m, the population values follow the trends of each whole population along the water column, with omnivores dominating at value 28, followed by herbivores at 3

Inferences from the parameters space a − z
For the parameters a (damage suffered by omnivores in their competition with herbivores) and z (attack rate by carnivores on omnivores) in all the seasons the herbivores, carnivores and omnivores show similar a trend. Along the water column the three populations decrease gradually although in different ways. Herbivores are dominant in the water layer up to depth 25m and decrease more in the deeper depths, while carnivores are substantially unchanged along all water column. Omnivores decrease, but in the deepest sites attain similar values as for the herbivores. These features are mainly influenced by two conditions: the reduction and absorption of the light, expressed by the model parameter ν, which increases with depth and the availability of food resources (phytoplankton for the herbivores), which is mainly responsible for the herbivores reproduction, the parameter r in (1). The ecosystem populations decrease with increasing depth occurs in different ways, depending on the season. During the winter in the deep stations at 50m and 75m, the herbivores decrease and the omnivores exceed the herbivores. This situation is repeated in the spring although on the surface the number of herbivores is quantitatively greater, due to larger food availability and primarily to the higher duration of the daily photoperiod, which supports the phytoplankton growth. Only at 75m below the surface, where the light absorption is maximal, the omnivores exceed herbivores. Carnivores are unchanged. Also in the summer, when the photoperiod is longer, the herbivores show the same features that they exhibit in the spring. During the fall instead, the population behaviors are similar to those of the winter.

Inferences from the parameters space j − z
Also for the parameters j (the losses of the omnivores population due to competition with carnivores) and z the factors most influencing the performance of the three populations of herbivores, omnivorous and carnivores appear to be the light and the availability of trophic resources.
In the shallow stations the omnivores are the predominant population. During the winter, in the deepest stations the carnivores decrease less than the omnivores and they are the predominant population, due to the reduction of the light and to the attack rate, that limit the presence of the omnivores. A similar situation occurs in the fall.
In the spring and above all during the summer, in the water column the temperature increases and this leads to thermal layers in the water body. The presence of thermoclines influences the phytoplankton distribution and consequently the herbivores distribution in the upper layers of the ocean up to 15m of depth, in which they are predominant. In the deepest station at 75m the herbivores decrease is caused by several concomitant factors, the decrease of available light, the presence of a colder thermocline and the reduction of phytoplankton. As a consequence, the omnivores become the predominant population.

Inferences from the parameters space
In the parameters b (attack rate by carnivores on herbivores), and g (damage suffered by carnivores in their competition with omnivores), the ecological conditions are similar with the other pairs of parameters discussed above. In fact the reduction and absorption of light increase with the depth, influencing the competition between omnivores and carnivores when the herbivores decrease. The herbivores population decreases mostly in the deepest two stations at 50m and 75m in all seasons, when the trophic resources are insufficient; this condition, together with the attack rate by carnivores on herbivores, favors the dominance of the omnivores mainly in the spring and summer.
During the cold seasons, in the fall and in the winter, the conditions in shallow layers are similar. However, they change dramatically when the trophic resources decrease, up to the point where the herbivores become close to vanishing, and a considerable reduction of the omnivores population is observed as the result of the competition between carnivores and omnivores.

Comments on the system behavior in some other parameter spaces
We now discuss some more situations in some of the remaining parameter spaces, not having the pretension of being exhaustive, in view of the too many possible pairwise combinations that arise in model (1), but above all because as already mentioned, the simulations relate on a particular, arbitrarily chosen, instant in time, so that they cannot be taken as a sure trend of the ecosystem, but only to illustrate some possible insights. In Figure 5 we report the three populations values in the spring for four of the five cases that we discuss below.
We simulated the populations behavior by combining the reciprocal competitions of herbivores and omnivores, namely the parameters a − c. Comparing the whole populations in the four seasons, an increase of the parameter a, i.e. the damage suffered by herbivores from omnivores, corresponds to an increase of the latter and a corresponding decrease of the former, as well as of the carnivores. This is observed in every season, although as discussed more thoroughly in the previous Section, the population values change at the various times of the year, as well as at the various depths in the water column. The parameter c, i.e. the damage suffered by omnivores, appears to affect less both herbivores and omnivores, while the carnivores appear to have their peak along an almost straight line, perhaps the bisectrix, of the parameter space.
Combined with the parameter g, i.e. damage incurred by carnivores in the competition with omnivores, instead, the parameter a shows a similar behavior, when at higher levels it enhances omnivores and hinders the herbivores. This behavior is qualitatively independent of the season, and not influenced by the parameter g. Carnivores experience a peak for intermediate values of a, which is decreasing for higher levels of the competition they suffer from omnivores. This very same behavior is found in the parameter space a − j, with the difference that the values of the three populations are not altered by the amount of losses incurred by omnivores in their competition with carnivores, j.
We now discuss the pair b − c, the direct hunting rate of carnivores on herbivores versus the omnivores losses due to competition with herbivores. In the range explored, the parameter b does not influence neither the herbivores nor the omnivores populations, while low levels of omnivores losses help this population to thrive and depress the herbivores. Carnivores too are hindered by a small value of c, they have a peak along a straight line almost parallel to the bisectrix of this parameter space, the maximum value increasing further away from the origin, i.e. with both b and c increasing.
We finally compare the hunting of carnivores on herbivores b with the damage inflicted on omnivores by the carnivores, j. Herbivores and omnivores have again opposite behavior, in this case however with high omnivores and low herbivores values for larger b and smaller j, while the reverse behavior occurs for smaller b and higher j. In each frame, left to right, are respectively shown the herbivores, omnivores and carnivores populations as functions of the given parameters.

Conclusion
In this paper, a mathematical model that describes the zooplankton spatio-temporal dynamics is analysed. The phytoplanktonic resource, not explicitly expressed in the model, is assumed to be unlimited. In the model, the influence of phytoplankton is expressed by the light absorption, along the entire water column. The light absorption values are obtained experimentally. The chemical-physical factors are neglected, therefore the phytoplankton dependence represents the main cause of spatial zooplanktonic populations dispersion. The model exhibits four feasible equilibrium points. For these points the stability conditions and the possible local bifurcations are determined. Moreover, from the eigenvalues study the point E 0 , that represents the absence of the three populations, i.e. ecosystem collapse, and the point E 4 , where the carnivorous zooplankton population is not present, are unstable. The points E 3 , E 5 and E 6 are stable if the conditions are verified (2), (7) and (10). For the analytical study of the local bifurcations, Sotomayor's theorem is used. From it, taking r as bifurcation parameter, in a neighborhood of the point E 3  From our numerical simulations, it is found that in the a − z parameter space, the three zooplankton popolations coexist at different levels, the herbivores dominating at every depth, followed by omnivores while carnivores are always found at low values. As expected, these populations decrease with the depth, consistently with the hypothesis that the reduced illumination hinders clorophyll production and therefore the phytoplankton abundance. In turn this affects primarily the herbivore population, that does not find nourishment, and the omnivores that feed also possibly on phytoplankton. As a consequence also the carnivores are affected, because the other reduced zooplankton populations cannot sustain them.
giving the explicit feasibility conditions r exp(−xν) > ℓ, s exp(−xν) + p > q (18) In case omnivores are absent, point E 5 , from the first equilibrium equation, solving for C we obtain a straight line that lies in the first quadrant only A > 0, if the left inequality in (3) is satisfied. From the third equilibrium equation we obtain directly This is a hyperbola raising up from the height at the origin kn −1 to the horizontal asymptote C = (k + eb)n −1 as H → +∞. A unique intersection among Φ and Ψ in the first quadrant therefore always exists if condition (3) holds.
These functions never intersect in the first quadrant if instead (4) is satisfied, while in the intermediate case either two or no intersections exist. In particular, equating (19) and (20), the quadratic in H is obtained: that has two roots if conditions (5) are verified. The population levels at this equilibrium are obtained from the roots of (21) and Φ(H), producing (6). For E 6 = (0, O 6 , C 6 ) from the second equilibrium equation, solving for C, we obtain while from the third one, once again solving for C, we find Note that χ(0) = B ∈ R and the function χ lies in the first quadrant only for M > 0, the right inequality in (8 (9), or none, in the opposite case.
At E 5± = (H 5 , 0, C 5± ) the Jacobian matrix factorizes to produce one explicit eigenvalue giving the first stability condition (7) while for the remaining minor J 5 using the first and third equilibrium equations we must have so that the Routh-Hurwitz conditions become the last two conditions in (7). In the case of the herbivore-free equilibrium E 6 one eigenvalue factors out, to provide the first stability condition in (10). Note further, that using the equilibrium equations, the diagonal entries in the Jacobian simplify as J 22 (E 6 ) = wz(1 + wO 6 ) −2 C 6 O 6 and J 33 (E 6 ) = −nC 6 , so that the Routh-Hurwitz conditions on the remaining minor provide the last two additional stability conditions in (10).

Bifurcations
Let f denote the right hand side of (1) and let ϕ = (H, O, C) T . For any nonzero vector v = (v 1 , v 2 , v 3 ) T and from the the Jacobian matrix Df = J, (24), we can evaluate Df T v as well as the following expressions and Near E 3 the following results hold. Theorem Using µ = r as bifurcation parameter, when it crosses the value µ 0 = r 0 in (11), near the equilibrium E 3 the model (1) has 1) no saddle-node bifurcation 2) a transcritical bifurcation 3) no pitchfork bifurcation.
Taking a different bifurcation parameter, we find the following result Theorem When the bifurcation parameter λ = k crosses the value λ 0 = k 0 of (12), near E 3 the model (1) shows 1) a saddle-node bifurcation 2) no transcritical bifurcation 3) no pitchfork bifurcation.
The local bifurcation analysis near E 5+ is summarized in the following result. Theorem When µ = r passes through the value µ 0 = r 0 of (13) near the equilibrium point E 5+ the model (1) has 1) no saddle-node bifurcation 2) no transcritical bifurcation 3) no pitchfork bifurcation.

Proof
The Jacobian (24) and its transpose evaluated at E 5+ and µ 0 have respectively the eigenvectors v = (1, 0, v 3  Since the conditions of Thm. 2.23 of [29] do not all hold, there are no local bifurcations at E 5 .