Skip to main content
Log in

Phytoplankton Competition for Nutrients and Light in a Stratified Lake: A Mathematical Model Connecting Epilimnion and Hypolimnion

  • Published:
Journal of Nonlinear Science Aims and scope Submit manuscript

Abstract

A mathematical model connecting epilimnion and hypolimnion is proposed to describe the competition of phytoplankton for nutrients and light in a stratified lake. The existence and stability of nonnegative steady-state solutions are completely characterized for all possible parameter ranges by means of stability analysis, bifurcation theory, and extensive simulations. The critical thresholds for settling speed of phytoplankton cells in the thermocline and the loss rate of phytoplankton are established, which determine the survival or extirpation of phytoplankton in epilimnion and hypolimnion. In particular, it is shown that in two extreme cases, the principle of competitive exclusion always holds in a stratified lake. We also consider the influence of environmental parameters on the vertical distribution and biomass density of phytoplankton via a systematic sensitivity analysis, and investigate their roles in phytoplankton blooms. These results can be used for the prediction of phytoplankton competition and blooms in a stratified lake.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9
Fig. 10
Fig. 11
Fig. 12
Fig. 13
Fig. 14
Fig. 15
Fig. 16

Similar content being viewed by others

References

  • Alijani, M.K., Wang, H., Elser, J.J.: Modeling the bacterial contribution to planktonic community respiration in the regulation of solar energy and nutrient availability. Ecol. Complex. 23, 25–33 (2015)

    Article  Google Scholar 

  • Boehrer, B., Schultze, M.: Stratification of lakes. Rev. Geophys. 46, 2 (2008)

    Article  Google Scholar 

  • Crandall, M.G., Rabinowitz, P.H.: Bifurcation from simple eigenvalues. J. Funct. Anal. 8(2), 321–340 (1971)

    Article  MathSciNet  MATH  Google Scholar 

  • Cusseddu, D., Edelstein-Keshet, L., Mackenzie, J.A., Portet, S., Madzvamuse, A.: A coupled bulk-surface model for cell polarisation. J. Theor. Biol. 481, 119–135 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Du, Y.H., Hsu, S.B.: Concentration phenomena in a nonlocal quasi-linear problem modelling phytoplankton. I. Existence. SIAM J. Math. Anal. 40(4), 1419–1440 (2008)

    Article  MathSciNet  MATH  Google Scholar 

  • Du, Y.H., Hsu, S.B.: Concentration phenomena in a nonlocal quasi-linear problem modelling phytoplankton. II. Limiting profile. SIAM J. Math. Anal. 40(4), 1441–1470 (2008)

    Article  MathSciNet  MATH  Google Scholar 

  • Du, Y.H., Hsu, S.B.: On a nonlocal reaction–diffusion problem arising from the modeling of phytoplankton growth. SIAM J. Math. Anal. 42(3), 1305–1333 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  • Edwards, A.M., Brindley, J.: Zooplankton mortality and the dynamical behaviour of plankton population models. Bull. Math. Biol. 61(2), 303–339 (1999)

    Article  MATH  Google Scholar 

  • Etemad-Shahidi, A., Imberger, J.: Anatomy of turbulence in thermally stratified lakes. Limnol. Oceanogr. 46(5), 1158–1170 (2001)

    Article  Google Scholar 

  • Gomez, D., Ward, M.J., Wei, J.C.: The linear stability of symmetric spike patterns for a bulk-membrane coupled Gierer–Meinhardt model. SIAM J. Appl. Dyn. Syst. 18(2), 729–768 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Grover, J.P.: Sink or swim? Vertical movement and nutrient storage in phytoplankton. J. Theor. Biol. 432, 38–48 (2017)

    Article  MATH  Google Scholar 

  • Heggerud, C.M., Wang, H., Lewis, M.A.: Transient dynamics of a stoichiometric cyanobacteria model via multiple-scale analysis. SIAM J. Appl. Math. 80(3), 1223–1246 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  • Henry, D.: Geometric theory of semilinear parabolic equations. Lecture Notes in Mathematics, vol. 840. Springer-Verlag, Berlin-New York (1981)

  • Hsu, S.B., Lou, Y.: Single phytoplankton species growth with light and advection in a water column. SIAM J. Appl. Math. 70(8), 2942–2974 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  • Hsu, S.B., Wang, F.B., Zhao, X.Q.: Global dynamics of zooplankton and harmful algae in flowing habitats. J. Differ. Equ. 255(3), 265–297 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  • Huisman, J., Weissing, F.J.: Light-limited growth and competition for light in well-mixed aquatic environments: an elementary mode. Ecology 75(2), 507–520 (1994)

    Article  Google Scholar 

  • Huisman, J., Arrayás, M., Ebert, U., Sommeijer, B.: How do sinking phytoplankton species manage to persist? Am. Nat. 159(3), 245–254 (2002)

    Article  Google Scholar 

  • Huisman, J., Arrayás, M., Ebert, U., Sommeijer, B.: Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum. Nature 439(7074), 322 (2006)

    Article  Google Scholar 

  • Jäger, C.G., Diehl, S.: Resource competition across habitat boundaries: asymmetric interactions between benthic and pelagic producers. Ecol. Monogr. 84(2), 287–302 (2014)

    Article  Google Scholar 

  • Jäger, C.G., Diehl, S., Emans, M.: Physical determinants of phytoplankton production, algal stoichiometry, and vertical nutrient fluxes. Am. Nat. 175(4), 91–104 (2010)

    Article  Google Scholar 

  • Jiang, J., Shen, A., Wang, H., Yuan, S.L.: Regulation of phosphate uptake kinetics in the bloom-forming dinoflagellates prorocentrum donghaiense with emphasis on two-stage dynamic process. J. Theor. Biol. 463, 12–21 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Jiang, D.H., Lam, K.Y., Lou, Y., Wang, Z.C.: Monotonicity and global dynamics of a nonlocal two-species phytoplankton model. SIAM J. Appl. Math. 79(2), 716–742 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Klausmeier, C.A., Litchman, E.: Algal games: the vertical distribution of phytoplankton in poorly mixed water columns. Limnol. Oceanogr. 46(8), 1998–2007 (2001)

    Article  Google Scholar 

  • Liu, P., Shi, J.P., Wang, Y.W.: Imperfect transcritical and pitchfork bifurcations. J. Funct. Anal. 251(2), 573–600 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  • Loladze, I., Kuang, Y., Elser, J.J.: Stoichiometry in producer-grazer systems: linking energy flow with element cycling. Bull. Math. Biol. 62(6), 1137–1162 (2000)

    Article  MATH  Google Scholar 

  • Lv, D.Y., Fan, M., Kang, Y., Blanco, K.: Modeling refuge effect of submerged macrophytes in lake system. Bull. Math. Biol. 78(4), 662–694 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  • Mei, L.F., Zhang, X.Y.: On a nonlocal reaction–diffusion–advection system modeling phytoplankton growth with light and nutrients. Discrete Contin. Dyn. Syst. Ser. B 17(1), 221–243 (2012)

    MathSciNet  MATH  Google Scholar 

  • Mei, L.F., Zhang, X.Y.: Existence and nonexistence of positive steady states in multi-species phytoplankton dynamics. J. Differ. Equ. 253(7), 2025–2063 (2012)

    Article  MathSciNet  MATH  Google Scholar 

  • Mischaikow, K., Smith, H., Thieme, H.R.: Asymptotically autonomous semiflows: chain recurrence and Lyapunov functions. Trans. Am. Math. Soc. 347(5), 1669–1685 (1995)

    Article  MathSciNet  MATH  Google Scholar 

  • Nie, H., Hsu, S.B., Wang, F.B.: Steady-state solutions of a reaction–diffusion system arising from intraguild predation and internal storage. J. Differ. Equ. 266(12), 8459–8491 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Nishimura, Y., Kim, C., Nagata, T.: Vertical and seasonal variations of bacterioplankton subgroups with different nucleic acid contents: possible regulation by phosphorus. Appl. Environ. Microbiol. 71(10), 5828–5836 (2005)

    Article  Google Scholar 

  • Pao, C.V.: Nonlinear Parabolic and Elliptic Equations. Plenum Press, New York (1992)

    MATH  Google Scholar 

  • Paquin-Lefebvre, F., Xu, B., DiPietro, K.L., Lindsay, A.E., Jilkine, A.: Pattern formation in a coupled membrane-bulk reaction–diffusion model for intracellular polarization and oscillations. J. Theor. Biol. 497(110242), 23 (2020)

    MathSciNet  Google Scholar 

  • Paquin-Lefebvre, F., Nagata, W., Ward, M.J.: Weakly nonlinear theory for oscillatory dynamics in a one-dimensional PDE-ODE model of membrane dynamics coupled by a bulk diffusion field. SIAM J. Appl. Math. 80(3), 1520–1545 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  • Peng, R., Zhao, X.Q.: A nonlocal and periodic reaction–diffusion–advection model of a single phytoplankton species. J. Math. Biol. 72(3), 755–791 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  • Ryabov, A.B., Rudolf, L., Blasius, B.: Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer. J. Theor. Biol. 263(1), 120–133 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  • Shi, J.P., Wang, X.F.: On global bifurcation for quasilinear elliptic systems on bounded domains. J. Differ. Equ. 246(7), 2788–2812 (2009)

    Article  MathSciNet  MATH  Google Scholar 

  • Song, D., Fan, M., Chen, M., Wang, H.: Dynamics of a periodic stoichiometric model with application in predicting and controlling algal bloom in Bohai Sea off China. Math. Biosci. Eng. 16(1), 119–138 (2019)

    Article  MathSciNet  Google Scholar 

  • Vasconcelos, F.R., Diehl, S., Rodríguez, P., Hedström, P., Karlsson, J., Byström, P.: Asymmetrical competition between aquatic primary producers in a warmer and browner world. Ecology 97(10), 2580–2592 (2016)

    Article  Google Scholar 

  • Wang, H., Smith, H.L., Kuang, Y., Elser, J.J.: Dynamics of stoichiometric bacteria–algae interactions in the epilimnion. SIAM J. Appl. Math. 68(2), 503–522 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  • Wang, F.B., Hsu, S.B., Zhao, X.Q.: A reaction–diffusion–advection model of harmful algae growth with toxin degradation. J. Differ. Equ. 259(7), 3178–3201 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  • Wang, Y., Shi, J.P., Wang, J.F.: Persistence and extinction of population in reaction–diffusion-advection model with strong allee effect growth. J. Math. Biol. 78(7), 2093–2140 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  • Wüest, A., Lorke, A.: Small-scale hydrodynamics in lakes. Annu. Rev. Fluid Mech. 35(1), 373–412 (2003)

    Article  MATH  Google Scholar 

  • Yoshiyama, K., Nakajima, H.: Catastrophic transition in vertical distributions of phytoplankton: alternative equilibria in a water column. J. Theor. Biol. 216(4), 397–408 (2002)

    Article  MathSciNet  Google Scholar 

  • Yoshiyama, K., Mellard, J.P., Litchman, E., Klausmeier, C.A.: Phytoplankton competition for nutrients and light in a stratified water column. Am. Nat. 174(2), 190–203 (2009)

    Article  MATH  Google Scholar 

  • Zagaris, A., Doelman, A.: Emergence of steady and oscillatory localized structures in a phytoplankton-nutrient model. Nonlinearity 24(12), 3437–3486 (2011)

    Article  MathSciNet  MATH  Google Scholar 

  • Zhang, J.M., Shi, J.P., Chang, X.Y.: A mathematical model of algae growth in a pelagic-benthic coupled shallow aquatic ecosystem. J. Math. Biol. 76(5), 1159–1193 (2018)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Junping Shi.

Additional information

Communicated by Mary Silber.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Partially supported by NSFC-11971088, NSFHLJ-LH2019A022, NSF-DMS-1853598, and NSERC Discovery Grant RGPIN-2020-03911 and Accelerator Grant RGPAS-2020-00090.

Appendices

Appendix A

In order to obtain the existence and stability of semi-trivial steady-state solutions of (2.1), we consider an eigenvalue problem

$$\begin{aligned} {\left\{ \begin{array}{ll} D\phi ''(x)=\lambda \phi (x),~x\in (0,x_h),\\ D\phi '(0)=\displaystyle \frac{\lambda }{\delta _1\lambda +\delta _2}\phi (0),~\phi (x_h)=0, \end{array}\right. } \end{aligned}$$
(A.1)

where \(D,\delta _1,\delta _2>0\). A direct calculation shows that if \(\lambda \) is a complex eigenvalue of (A.1), then \({\bar{\lambda }}\) is also a complex eigenvalue of (A.1).

Lemma A.1

If \(\lambda \) is an eigenvalue of (A.1), then \({{\,\mathrm{Re}\,}}\lambda <0\).

Proof

If \(\lambda \) is an eigenvalue (real or complex valued) and \(\phi (x)\) is the corresponding eigenfunction, then

$$\begin{aligned} \lambda \displaystyle \int _0^{x_h}|\phi (x)|^2\mathrm{d}x= D\displaystyle \int _0^{x_h}\phi ''(x){\bar{\phi }}(x)\mathrm{d}x =-\displaystyle \frac{\lambda |\phi (0)|^2}{\delta _1\lambda +\delta _2}-D\displaystyle \int _0^{x_h} |\phi '(x)|^2\mathrm{d}x. \end{aligned}$$

Let

$$\begin{aligned} c_1=\displaystyle \int _0^{x_h}|\phi (x)|^2\mathrm{d}x,~c_2=D\displaystyle \int _0^{x_h} |\phi '(x)|^2\mathrm{d}x. \end{aligned}$$

Then, \(\lambda \) satisfies

$$\begin{aligned} \delta _1c_1\lambda ^2+(\delta _2c_1+|\phi (0)|^2+\delta _1c_2)\lambda +\delta _2c_2=0. \end{aligned}$$

This means that \({{\,\mathrm{Re}\,}}\lambda <0\). This completes the proof. \(\square \)

Proof of Theorem 3.1

From (3.1), it is easy to see that \(E_1\equiv (0,M_b,0,M_b)\) is the unique nutrients-only semi-trivial steady state solution of (2.1). The stability of \(E_1\) is determined by the eigenvalue problem

$$\begin{aligned}&\lambda \xi =\left( \delta _a^*-\delta \right) \xi , \end{aligned}$$
(A.2a)
$$\begin{aligned}&\lambda \zeta =\theta \left( p\delta -\delta _a^*-\displaystyle \frac{a}{x_e}\right) \xi -\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0),\end{aligned}$$
(A.2b)
$$\begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+ \left( rf(M_b)g(I(x,0,0)) -\delta \right) \varphi ,~0<x<x_h,\end{aligned}$$
(A.2c)
$$\begin{aligned}&\lambda \psi =D_m\psi ''(x)+ \left( \theta p\delta -\theta rf(M_b)g(I(x,0,0))\right) \varphi ,~0<x<x_h,\end{aligned}$$
(A.2d)
$$\begin{aligned}&D_b\varphi '(0)-v\varphi (0)=-a\xi ,~D_b\varphi '(x_h)-v\varphi (x_h)=0, \end{aligned}$$
(A.2e)
$$\begin{aligned}&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0. \end{aligned}$$
(A.2f)

To establish the local stability of \(E_1\), we let \(\lambda _1\) be the eigenvalue of (A.2) with the largest real part and let \((\xi ,\zeta ,\varphi ,\psi )\) be the corresponding eigenfunction. Note that the linearized system (A.2) is partially decoupled. We consider the following three cases: (i) \(\xi \ne 0\); (ii) \(\xi =0, \; \varphi \not \equiv 0\); or (iii) \(\xi =0, \varphi \equiv 0\).

  1. Case (i):

    \(\xi \ne 0\). In this case, the stability of \(E_1\) is determined by (A.2a). Then, \(\lambda _1=\delta _a^*-\delta \).

  2. Case (ii):

    \(\xi =0, \; \varphi \not \equiv 0\). This case means that the stability of \(E_1\) is determined by (A.2c) and its boundary condition (A.2e). Then, \(\lambda _1=\delta _*-\delta \).

  3. Case (iii):

    \(\xi =0, \; \varphi \equiv 0\). In this case, (A.2) reduces to

    $$\begin{aligned} {\left\{ \begin{array}{ll} \lambda \zeta =-\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0),&{}\\ \lambda \psi (x)=D_m\psi ''(x),~~0<x<x_h,&{}\\ D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0.&{} \end{array}\right. } \end{aligned}$$
    (A.3)

    From the boundary condition of (A.3), we have \(\zeta =0\) if \(\psi \equiv 0\). On the other hand, if \(\zeta =0\) in (A.3), we also have \(\psi \equiv 0\). This shows that \(\zeta \ne 0\) and \(\psi (x)\not \equiv 0\). Then, either \(\lambda =-b/x_e=-D_mk^2\pi ^2/x_h^2\) for some \(k\in {\mathbb {N}}\) with the corresponding eigenfunction \((\zeta ,\psi )=(1,-\sqrt{bx_e/D_m}\sin (\sqrt{b/(x_eD_m)}x)\); or \(\lambda \ne -b/x_e\), \(\displaystyle \zeta =\displaystyle \frac{b}{x_e\lambda +b}\psi (0)\), and \(\psi \) satisfies

    $$\begin{aligned} {\left\{ \begin{array}{ll} \lambda \psi (x)=D_m\psi ''(x),~~0<x<x_h,&{}\\ D_m\psi '(0)=\displaystyle \frac{bx_e\lambda }{x_e\lambda +b}\psi (0),~\psi (x_h)=0.&{} \end{array}\right. } \end{aligned}$$
    (A.4)

It follows from Lemma A.1 that all eigenvalues of (A.4) have negative real parts. Therefore, we conclude that in case (iii), \({{\,\mathrm{Re}\,}}\lambda _1\) is negative.

Summarizing the above cases (i)-(iii), we conclude that if (3.6) holds, then \({{\,\mathrm{Re}\,}}\lambda _1<0\) and \(E_1\) is locally asymptotically stable; Conversely, if (3.7) holds, then \({{\,\mathrm{Re}\,}}\lambda _1>0\) and \(E_1\) is unstable.

Next we prove that if (3.8) holds, then \(E_1\) is globally asymptotically stable. From the first equation of (2.1), we have

$$\begin{aligned} \displaystyle \frac{dA(t)}{dt}\le rA\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x -\displaystyle \frac{a}{x_e}A-\delta A, \end{aligned}$$

which implies that \(\displaystyle \lim _{t\rightarrow \infty }A(t)=0\) if \(\delta >\delta _a^{**}\). From the theory of asymptotical autonomous systems (Mischaikow et al. 1995), the third equation of (2.1) reduces to a limiting system

$$\begin{aligned} \displaystyle \frac{\partial B(x,t)}{\partial t}&=D_b\displaystyle \frac{\partial ^2B}{\partial x^2}-v\displaystyle \frac{\partial B}{\partial x}+rBf(M)g(I(x,A,B))-\delta B\\&\le D_b\displaystyle \frac{\partial ^2B}{\partial x^2}-v\displaystyle \frac{\partial B}{\partial x}+rg(I(0,0,0))B-\delta B,~0<x<x_h,~t>0,\\ D_b\displaystyle \frac{\partial B(0,t)}{\partial x}&-vB(0,t)= D_b\displaystyle \frac{\partial B(x_h,t)}{\partial x}-vB(x_h,t)=0,~t>0. \end{aligned}$$

This implies that B(xt) converges to 0 uniformly for \(x\in [0,x_h]\) as \(t\rightarrow \infty \) if \(\delta >\delta _{**}\) by the comparison principle of parabolic equations. By applying the theory of asymptotical autonomous systems again, (2.1) reduces to a limiting system

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{dN(t)}{dt}=\displaystyle \frac{b}{x_e}(M(0,t)-N(t)),&{}t>0,\\ \displaystyle \frac{\partial M(x,t)}{\partial t}=D_m\displaystyle \frac{\partial ^2M}{\partial x^2},&{}0<x<x_h,~t>0,\\ D_m\displaystyle \frac{\partial M(0,t)}{\partial x}=b(M(0,t)-N(t)),~M(x_h,t)=M_b,&{}t>0. \end{array}\right. } \end{aligned}$$
(A.5)

Define Lyapunov functional \(V:{\mathbb {R}}\times C([0,L])\rightarrow {\mathbb {R}}\) by

$$\begin{aligned} V(N,M)=\displaystyle \frac{x_e}{2}(N-M_b)^2+\displaystyle \frac{1}{2}\displaystyle \int _0^{x_h}(M(z)-M_b)^2\mathrm{d}z. \end{aligned}$$

For an arbitrary solution (N(t), M(xt)) of (A.5) with nonnegative initial values, we obtain

$$\begin{aligned} \displaystyle \frac{dV(N(t),M(x,t))}{dt}&=x_e(N(t)-M_b)\displaystyle \frac{d N}{dt} +\displaystyle \int _0^{x_h}(M(z,t)-M_b)\displaystyle \frac{\partial M}{\partial t}\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))+D_m\displaystyle \int _0^{x_h}(M(z,t)-M_b)\displaystyle \frac{\partial ^2 M}{\partial z^2}\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))\\&\quad +D_m\displaystyle \frac{\partial M}{\partial z}(M(z,t)-M_b)\Big |_0^{x_h}-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\\&=b(N(t)-M_b)(M(0,t)-N(t))\\&\quad -b(M(0,t)-N(t))(M(0,t)-M_b)-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\\&=-b(N(t)-M(0,t))^2-\displaystyle \int _0^{x_h} \left( \displaystyle \frac{\partial M}{\partial z}\right) ^2\mathrm{d}z\le 0. \end{aligned}$$

Note that \(dV(\cdot )/dt=0\) holds if and only if \(\partial M/\partial z=0\) and \(N(t)=M(0,t)\). It follows from \(M(x_h,t)=M_b\) that \(N(t)\equiv M(x,t)\equiv M_b\). By using the LaSalle’s Invariance Principle (Henry 1981), we conclude that (N(t), M(xt)) converges to \((M_b,M_b)\) uniformly for \(x\in [0,x_h]\) as \(t\rightarrow \infty \), and \(E_1\) is globally asymptotically stable for (2.1) with respect to any nonnegative initial value. \(\square \)

Proof of Theorem 3.3

The steady-state equation (3.2) can be explicitly solved. The equation of N implies that \(b(M(0)-N)=(1-p)\theta \delta x_eA\). Combining the equation of M with its boundary conditions \(D_mM'(0)=(1-p)\theta \delta x_eA\) and \(M(x_h)=M_b\) gives

$$\begin{aligned} M(x)=M_b-\displaystyle \frac{x_h-x}{D_m}(1-p)\theta \delta x_eA, \end{aligned}$$

and then

$$\begin{aligned} N=M_b-\left( \displaystyle \frac{x_h}{D_m}+\displaystyle \frac{1}{b}\right) (1-p)\theta \delta x_eA. \end{aligned}$$

If \(p=1\) and (3.9) holds, then \(M(x)=N=M_b\), and there exists a unique positive \(A_2^*\) satisfying \(\displaystyle f(M_b)\displaystyle \frac{r}{x_e}\int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x=\delta \) since g is a strictly decreasing function of A. If \(p\in [0,1)\), then we let

$$\begin{aligned} Q=\delta A,~~\beta =\left( \displaystyle \frac{x_h}{D_m}+\displaystyle \frac{1}{b}\right) (1-p)\theta x_e. \end{aligned}$$

From the first equation of (3.2), we consider

$$\begin{aligned} u(\delta ,Q)= & {} \displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g\left( I\left( x,\displaystyle \frac{Q}{\delta }\right) \right) \mathrm{d}x -\delta \left( 1+\displaystyle \frac{\gamma _1}{M_b-\beta Q}\right) ~\text{ for }~(\delta ,Q)\in (0,\infty )\\&\quad \times \left( 0,\displaystyle \frac{M_b}{\beta }\right) . \end{aligned}$$

A direct calculation gives \(\partial u/\partial Q<0\). Note that \(\displaystyle \lim _{Q\rightarrow 0^+}u(\delta ,Q)=(1+\gamma _1/M_b)(\delta ^{*}_0-\delta )>0\) if \(0<\delta <\delta ^{*}_0\) and \(\displaystyle \lim _{Q\rightarrow (M_b/\beta )^{-}}u(\delta ,Q)=-\infty \). Then for any fixed \(\delta \in (0,\delta _0^*)\), there exists a unique positive \(Q_\delta \) such that \(u(\delta ,Q_\delta )=0\), and \(A_2^*=Q_\delta /\delta \) is desired unique solution. This proves part (i).

Next we establish the stability of \(E_2\). The stability of \(E_2\) is determined by the eigenvalue problem

$$\begin{aligned}&\lambda \xi =-h_1\xi +h_2\zeta , \end{aligned}$$
(A.6a)
$$\begin{aligned}&\lambda \zeta = (-(1-p)\theta \delta +\theta h_1)\xi -\left( \displaystyle \frac{b}{x_e}+\theta h_2\right) \zeta +\displaystyle \frac{b}{x_e}\psi (0),\end{aligned}$$
(A.6b)
$$\begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+ \left( h_3(x) -\delta \right) \varphi ,~0<x<x_h,\end{aligned}$$
(A.6c)
$$\begin{aligned}&\lambda \psi =D_m\psi ''(x)+ \left( p\theta \delta -\theta h_3(x)\right) \varphi ,~0<x<x_h,\end{aligned}$$
(A.6d)
$$\begin{aligned}&D_b\varphi '(0)-v\varphi (0)=D_b\varphi '(x_h)-v\varphi (x_h)=0,\end{aligned}$$
(A.6e)
$$\begin{aligned}&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0, \end{aligned}$$
(A.6f)

where

$$\begin{aligned} \begin{aligned}&h_1=\displaystyle \frac{r\gamma _2\mu (A^*_2)N_2^*A^*_2}{N_2^*+\gamma _1},~ h_2=\displaystyle \frac{r\gamma _1A^*_2}{(N_2^*+\gamma _1)^2}\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x,\\&h_3(x)=rf(M_2^*)g(I(x,A^*_2,0)) \end{aligned} \end{aligned}$$
(A.7)

and

$$\begin{aligned} \mu (A^*_2)=\displaystyle \frac{l\left( \ln \displaystyle \frac{I_{in}+\gamma _2}{I_{in}e^{-x_e(K_{bg}+lA_2^*)}+\gamma _2}-\displaystyle \frac{I_{in}x_e(K_{bg}+lA^*_2)e^{-x_e(K_{bg}+lA^*_2)}}{I_{in}e^{-x_e(K_{bg}+lA^*_2)}+\gamma _2} \right) }{x_e(K_{bg}+lA^*_2)^2}>0.\nonumber \\ \end{aligned}$$
(A.8)

Again let \(\lambda _1\) be the eigenvalue of (A.6) with largest real part, and let \((\xi ,\zeta ,\varphi ,\psi )\) be the corresponding eigenfunction. Note that the eigenvalue problem (A.6) is partially decoupled. We consider the following two cases: (i) \(\varphi \not \equiv 0\); (ii) \(\varphi \equiv 0\).

Case (i): \(\varphi \not \equiv 0\). Then, the stability of \(E_2\) is determined by the subsystem

$$\begin{aligned} \begin{aligned}&\lambda \varphi =D_b\varphi ''(x)-v\varphi '(x)+\left( h_3(x)-\delta \right) \varphi ,~0<x<x_h,\\&D_b\varphi '(0)-v\varphi (0)=D_b\varphi '(x_h)-v\varphi (x_h)=0, \end{aligned} \end{aligned}$$
(A.9)

with \(\xi ,\zeta ,\psi \) solved from (A.6a), (A.6b), (A.6d) and (A.6f). Then, \(\lambda _1=\lambda _1(D_b,v,h_3(x)-\delta )=\lambda _1(D_b,v,h_3(x))-\delta \).

Case (ii): \(\varphi \equiv 0\). Then, \((\xi ,\zeta ,\psi )\) satisfies

$$\begin{aligned} \begin{aligned}&\lambda \xi =-h_1\xi +h_2\zeta ,\\&\lambda \zeta = (-(1-p)\theta \delta +\theta h_1)\xi -\left( \displaystyle \frac{b}{x_e}+\theta h_2\right) \zeta +\displaystyle \frac{b}{x_e}\psi (0),\\&\lambda \psi =D_m\psi ''(x),~0<x<x_h,\\&D_m\psi '(0)=b(\psi (0)-\zeta ),~\psi (x_h)=0. \end{aligned} \end{aligned}$$
(A.10)

By using similar arguments as those in Theorem 3.1 and Lemma A.1, we conclude that all eigenvalues of (A.10) have negative real parts. Based on the analysis above, \(E_2\) is locally asymptotically stable with respect to (2.1) if (3.11) holds, while \(E_2\) is unstable if (3.12) holds.

It can be shown that \(A_2^*\) is strictly decreasing in \(\delta \) and \(\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}A_2^*=0\). This also implies that \(\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}N_2^*=M_b\) and \(\displaystyle \lim _{\delta \rightarrow \delta _0^{*-}}M_2^*(x)=M_b\). Hence, when \(\delta \rightarrow \delta _0^{*-}\), we have

$$\begin{aligned} \begin{aligned}&\lim _{\delta \rightarrow \delta _0^{*-}}\left[ \lambda _1\left( D_b,v,rf(M_2^*)g(I(x,A_2^*,0) \right) -\delta \right] \\ =&\lim _{\delta \rightarrow \delta _0^{*-}}\left[ \lambda _1\left( D_b,v,rf(M_2^*)g(I(x,A_2^*,0) \right) \right. \\&\left. -\lambda _1\left( D_b,v,f(N_2^*)\displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x\right) \right] \\ =&\lambda _1\left( D_b,v,rf(M_b)g(I(x,0,0)\right) -\lambda _1\left( D_b,v,f(M_b)\displaystyle \frac{r}{x_e} \displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x\right) <0. \end{aligned} \end{aligned}$$

Therefore, (3.11) holds and \(E_2\) is locally asymptotically stable when \(\delta \) is sufficiently close to \(\delta _0^*\). When \(p=1\), \(E_2\equiv (A_2^*,M_b,0,M_b)\) where \(A_2^*\) satisfies (3.10) with \(N_2^*=M_b\). From (3.5) and the monotonicity of the principal eigenvalue on the weight functions, we have

$$\begin{aligned} \lambda _1\left( D_b,v,rf(M_b)g(I(x,A_2^*,0)\right) < \lambda _1\left( D_b,v,f(M_b)\displaystyle \frac{r}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_2^*))\mathrm{d}x\right) =\delta , \end{aligned}$$

so (3.11) holds and \(E_2\) is locally asymptotically stable for any \(0<\delta <\delta _0^*\). \(\square \)

To obtain the global existence of \(E_3\), we first establish the following a priori estimates for nonnegative solutions \((N_3,B_3(x),M_3(x))\) of (3.3).

Lemma A.2

Assume that \((N_3,B_3(x),M_3(x))\in {\mathbb {R}}^+\times Y^2\) is a nonnegative solution of (3.3) with \(N_3>0,B_3,M_3\not \equiv 0\). Then,

  1. (i)

    \(0<\delta <\delta _{**}\), where \(\delta _{**}\) is defined in (3.5);

  2. (ii)

    For any \(\epsilon >0\), there exists a positive constant \(K(\epsilon )\) such that \(0<B_3(x)\le K(\epsilon )\) and \(0<N_3,M_3(x)\le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m\) on \([0,x_h]\) for \(\delta \in [\epsilon ,\delta _{**})\).

Proof

(i) Let \(R(x)=B_3(x)e^{-(v/D_b)x}\). Then, R(x) satisfies

$$\begin{aligned} {\left\{ \begin{array}{ll} -D_bR''(x)-vR'(x)+\delta R=rf(M_3)g(I(x,0,B_3))R\ge 0,&{}0<x<x_h,\\ R'(0)=0,~R'(x_h)=0.&{} \end{array}\right. } \end{aligned}$$

It follows from the strong maximum principle that \(R(x)>0\) on \([0,x_h]\), and consequently, \(B_3(x)>0\) on \([0,x_h]\). From (3.3), we have

$$\begin{aligned} \lambda _1\left( D_b,v,rf(M_3)g(I(x,0,B_3))\right) =\delta \end{aligned}$$

with corresponding principal eigenfunction \(B_3\). It follows from the monotonicity of principal eigenvalue with respect to the weight functions that

$$\begin{aligned} \delta =\lambda _1\left( D_b,v,rf(M_3)g(I(x,0,B_3))\right) <\lambda _1\left( D_b,v,rg(I(0,0,0))\right) =\delta _{**}. \end{aligned}$$

(ii) Fix \(\epsilon >0\). If \(B_3\) is not bounded for all \(\delta \in [\epsilon ,\delta _{**})\), then there are a sequence \(\delta ^i\in [\epsilon ,\delta _{**})\) and corresponding positive solutions \((N_3^i,B_3^i(x),M^i_3(x))\) such that \(\Vert B_3^i\Vert _\infty \rightarrow \infty \) and \(\delta ^i\rightarrow \delta ^0\in [\epsilon ,\delta _{**}]\) as \(i\rightarrow \infty \). Let \(\mu _i=B_3^i/\Vert B_3^i\Vert _\infty \). From the equation of \(B_3\) in (3.3), we get

$$\begin{aligned} {\left\{ \begin{array}{ll} -D_b\mu _i''(x)+v\mu _i'(x)=[rf(M_3^i)g(I(x,0,B_3^i))-\delta ^i] \mu _i(x)=0,&{}0<x<x_h,\\ D_b\mu _i'(0)-v\mu _i(0)=D_b\mu _i'(x_h)-v\mu _i(x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$
(A.11)

Since the right hand side of (A.11) is uniformly bounded, by using \(L^p\) theory for elliptic operators and by passing to a subsequence, we obtain that \(\mu _i\rightarrow \mu \) in \(W^{2,p}([0,x_h)]\) (and also in \(C^{1,\alpha }([0,x_h])\) from Sobolev’s embedding) as \(i\rightarrow \infty \). Since \(rf(M_3^i)g(I(x,0,B_3^i))\) is bounded in \(L^{\infty }([0,x_h])\), we may assume that \(rf(M_3^i)g(I(x,0,B_3^i))\rightarrow l^0\) weakly in \(L^2([0,x_h])\). Hence, \(\mu \) satisfies (in the weak sense)

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\mu ''(x)-v\mu '(x)+(l^0-\delta ^0) \mu (x)=0,&{}0<x<x_h,\\ D_b\mu '(0)-v\mu (0)=D_b\mu '(x_h)-v\mu (x_h)=0.&{} \end{array}\right. } \end{aligned}$$
(A.12)

From the strong maximum principle, \(\mu (x)>0\) on \([0,x_h]\) since \(\mu \ge 0\) and \(\Vert \mu \Vert _\infty =1\). This implies that \(B_3^i=||B_3^i||_{\infty }\mu _i\rightarrow \infty \) uniformly on \([0,x_h]\), and thus, \(l^0=0\). Integrating (A.12) we obtain \(0=\delta ^0\displaystyle \int _0^{x_h}\mu (x)\mathrm{d}x>0.\) That is a contradiction. Therefore, there is a positive constant \(K(\epsilon )\) such that \(0<B_3(x)\le K(\epsilon )\) on \([0,x_h]\) for all \(\delta \in [\epsilon ,\delta _{**})\).

From the strong maximum principle, we have \(M_3(x)>0\) on \([0,x_h]\). For any \(x\in [0,x_h]\), we have

$$\begin{aligned} |D_mM'_3(x)|&=\left| D_m\displaystyle \int _0^{x}M''_3(s)ds\right| = \theta \left| \displaystyle \int _0^x\left( rf(M_3)g(I(x,0,B_3))-p\delta \right) B_3ds\right| \\&\le \theta x_h(r+p\delta _{**})K(\epsilon ) \end{aligned}$$

and consequently

$$\begin{aligned} \begin{aligned} |M_3(x)|&=|M_3(x_h)+M_3(x)-M_3(x_h)|\le |M_3(x_h)|+|M_3(x)-M_3(x_h)|\\&\le M_b+\left| \displaystyle \int _x^{x_h}M_3'(s)ds \right| \le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m. \end{aligned} \end{aligned}$$
(A.13)

From \(M_3(0)=N_3\), we have \(0<N_3\le M_b+\theta x_h^2(r+p\delta _{**})K(\epsilon )/D_m\). \(\square \)

Proof of Theorem 3.5

We prove the existence of \(E_3\) using bifurcation theory and show that the solution \(E_3\) bifurcates from the line of nutrient-only semi-trivial steady state \(E_1\) at \(\delta =\delta _{*}\) with parameter \(\delta \). We first use local bifurcation theory in Crandall and Rabinowitz (1971) to show that \(E_3\) bifurcates from the line of \(E_1\) at \(\delta =\delta _{*}\).

Recall function spaces \(X_1,X_2\) defined in (3.13) and define \(Y=C([0,x_h])\). We define a nonlinear mapping \(F:{\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2\rightarrow {\mathbb {R}}\times Y^2\times {\mathbb {R}}^2\) by

$$\begin{aligned} F(\delta ,N,B,M)=\begin{pmatrix} M(0)-N\\ D_bB''-vB'+rBf(M)g(I(x,0,B))-\delta B\\ D_mM''+\theta p\delta B-\theta rBf(M)g(I(x,0,B))\\ D_mM'(0)-b(M(0)-N)\\ M(x_h)-M_b \end{pmatrix}. \end{aligned}$$

It is easy to see that \(F(\delta ,M_b,0,M_b)=0\). Let \(H:=F_{(N,B,M)}(\delta _*,M_b,0,M_b)\) be the Frechét derivative of F with respect to (NBM) at \((\delta _*,M_b,0,M_b)\). For any \((\zeta ,\varphi ,\psi )\in {\mathbb {R}}\times X_1\times X_2\), we have

$$\begin{aligned} \begin{aligned} H[\zeta ,\varphi ,\psi ] =\begin{pmatrix}\psi (0)-\zeta \\ D_b\varphi ''(x)-v\varphi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi \\ D_m\psi ''(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi \\ D_m\psi '(0)-b(\psi (0)-\zeta )\\ \psi (x_h) \end{pmatrix}. \end{aligned} \end{aligned}$$
(A.14)

If \((\zeta _1,\varphi _1,\psi _1)\in \ker H\), then

$$\begin{aligned}&\psi _1(0)-\zeta _1=0, \end{aligned}$$
(A.15a)
$$\begin{aligned}&D_b\varphi ''_1(x)-v\varphi '_1(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi _1=0,~0<x<x_h, \end{aligned}$$
(A.15b)
$$\begin{aligned}&D_m\psi ''_1(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi _1=0,~0<x<x_h, \end{aligned}$$
(A.15c)
$$\begin{aligned}&\psi '_1(0)=0,~\psi _1(x_h)=0. \end{aligned}$$
(A.15d)

Recall that \(\delta _*\) is the principal eigenvalue of (3.4) with \(q(x)=rf(M_b)g(I(x,0,0))\), and let \({\bar{\varphi }}\in X_1\) be the corresponding positive eigenfunction for the principal eigenvalue \(\delta _*\). Then, \({\bar{\varphi }}\) is the unique solution of (A.15b) up to a constant multiple. And there exist unique functions \({\bar{\psi }}\in X_2\) and \({\bar{\zeta }}\in {\mathbb {R}}\) satisfying (A.15c), (A.15d) and (A.15a). Hence, \(\dim \ker H=1\) and \(\ker H={{\,\mathrm{span}\,}}\{({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})\}.\)

If \((\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\in {{\,\mathrm{range}\,}}H\), then there exists \((\zeta _2,\varphi _2,\psi _2)\in {\mathbb {R}}\times X_1\times X_2\) such that

$$\begin{aligned} \begin{aligned}&\psi _2(0)-\zeta _2=\sigma _1,\\&D_b\varphi ''_2(x)-v\varphi '_2(x)+\left( rf(M_b)g(I(x,0,0))-\delta _*\right) \varphi _2=\sigma _2(x),~0<x<x_h,\\&D_m\psi ''_2(x)+\left( \theta p\delta _*-\theta rf(M_b)g(I(x,0,0))\right) \varphi _2=\sigma _3(x),~0<x<x_h,\\&D_m\psi '_2(0)-b(\psi _2(0)-\zeta _2)=\sigma _4,~\psi _2(x_h)=\sigma _5. \end{aligned} \end{aligned}$$
(A.16)

Multiplying both sides of (A.15b) and the second equation of (A.16) by \(\varphi _2e^{-(v/D_b)x}\) and \(\varphi _1e^{-(v/D_b)x}\), respectively, subtracting, and integrating on \([0,x_h]\), we obtain

$$\begin{aligned}&\displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}\varphi _1(x)\mathrm{d}x\\=&D_b\displaystyle \int _0^{x_h}\left( \left( \varphi '_2(x) e^{-(v/D_b)x}\right) '\varphi _1(x)-\left( \varphi '_1(x) e^{-(v/D_b)x}\right) '\varphi _2(x)\right) \mathrm{d}x\\ =&D_b\varphi '_2(x)e^{-(v/D_b)x}\varphi _1(x)\Big |_0^{x_h}-D_b\varphi '_1(x)e^{-(v/D_b)x} \varphi _2(x)\Big |_0^{x_h}=0. \end{aligned}$$

This implies that

$$\begin{aligned} \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x=0 \end{aligned}$$

and

$$\begin{aligned} {{\,\mathrm{range}\,}}H=\left\{ (\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\in {\mathbb {R}}\times Y^2\times {\mathbb {R}}^2:\displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x=0\right\} . \end{aligned}$$

Then, \({{\,\mathrm{codim}\,}}{{\,\mathrm{range}\,}}H=1\). Moreover, we also have

$$\begin{aligned} F_{\delta (N,B,M)}(\delta _*,M_b,0,M_b) ({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})=(0,-{\bar{\varphi }},p\theta {\bar{\varphi }},0,0), \end{aligned}$$

so \(F_{\delta ,(N,B,M)}(\delta _*,M_b,0,M_b)({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}) \not \in {{\,\mathrm{range}\,}}H\) as \(\displaystyle \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}^2(x)\mathrm{d}x\ne 0\).

From Theorem 1.7 in Crandall and Rabinowitz (1971), there exists a positive constant \(\varepsilon _3>0\) such that all solutions of (3.3) near \((\delta _*,M_b,0,M_b)\) lie on a smooth curve

$$\begin{aligned} \Gamma _3=\{(\delta _3(s),N_3^*(s),B_3^*(s,x),M_3^*(s,x)):0<s<\varepsilon _3\} \end{aligned}$$

with the form

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}N_3^*(s)=M_b+s{\bar{\zeta }}+o(s),\;\; B_3^*(s,x)=s{\bar{\varphi }}(x)+o(s),\\ &{}M_3^*(s,x)=M_b+s{\bar{\psi }}(x)+o(s). \end{array}\right. } \end{aligned}$$
(A.17)

Let \({\bar{l}}\) be a linear functional on \({\mathbb {R}}\times Y^2\times {\mathbb {R}}^2\) by

$$\begin{aligned} \langle {\bar{l}},(\sigma _1,\sigma _2,\sigma _3,\sigma _4,\sigma _5)\rangle = \displaystyle \int _0^{x_h}\sigma _2(x)e^{-(v/D_b)x}{\bar{\varphi }}(x)\mathrm{d}x. \end{aligned}$$

Then, from Liu et al. (2007), we have

$$\begin{aligned} \begin{aligned} \delta _3'(0)&=-\displaystyle \frac{\left\langle {\bar{l}},F_{(N,B,M)(N,B,M)}\left( \delta _*, M_b,0,M_b\right) [{\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}]^2\right\rangle }{2\left\langle {\bar{l}}, F_{\delta ,(N,B,M)}\left( \delta _*, M_b,0,M_b\right) [{\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }}]\right\rangle }\\&=-\displaystyle \frac{\displaystyle \int _0^{x_h}\displaystyle \frac{rI(x,0,0)e^{-(v/D_b)x}}{(M_b+\gamma _1)(I(x,0,0)+\gamma _2)} \mu (x,{\bar{\varphi }},{\bar{\psi }})\mathrm{d}x}{\displaystyle \int _0^{x_h}e^{-(v/D_b)x} {\bar{\varphi }}^2(x)\mathrm{d}x}, \end{aligned} \end{aligned}$$
(A.18)

where

$$\begin{aligned} \mu (x,{\bar{\varphi }},{\bar{\psi }})= \displaystyle \frac{\gamma _2 lM_b}{I(x,0,0)+\gamma _2}{\bar{\varphi }}^2(x)\displaystyle \int _0^x {\bar{\varphi }}(s)ds- \displaystyle \frac{\gamma _1}{M_b+\gamma _1}{\bar{\varphi }}^2(x){\bar{\psi }}(x),~x\in [0,x_h]. \end{aligned}$$

We claim that \({\bar{\psi }}(x)\le 0\) on \(x\in [0,x_h]\). Let \(y(x)=rf(M_b)g(I(x,0,0)).\) Then y(x) is a strictly decreasing function on \(x\in [0,x_h]\). It follows from (A.15b)-(A.15d) that \({\bar{\psi }}'(x_h)\ge 0\). From (A.15c), we have \({\bar{\psi }}''(x)>0\) on \(x\in [0,x_h]\) if \(p\le y(x_h)/\delta _*\). Combining with its boundary conditions \({\bar{\psi }}'(x)=0,{\bar{\psi }}(x_h)=0\), we get \({\bar{\psi }}(x)\le 0\) on \(x\in [0,x_h]\). If \(p>y(x_h)/\delta _*\), then there is a unique \(x_0\in (0,x_h)\) such that \({\bar{\psi }}''(x)>0\) on \(x\in (0,x_0)\) and \({\bar{\psi }}''(x)<0\) on \(x\in (x_0,x_h)\). This means that \({\bar{\psi }}'(x)\) is a strictly increasing function on \(x\in [0,x_0)\) and a strictly decreasing function on \(x\in (x_0,x_h]\). Hence \({\bar{\psi }}'(x)>0\) on \(x\in (0,x_h)\) since \({\bar{\psi }}'(0)=0\) and \({\bar{\psi }}'(x_h)\ge 0\). It follows from \({\bar{\psi }}(x_h)=0\) that \({\bar{\psi }}(x)\le 0\) on \(x\in [0,x_h]\). From (A.18), we have \(\delta _3'(0)<0\), which implies that the bifurcation at \((\delta _*,M_b,0,M_b)\) is backward. This completes the proof of part (iii).

Now we turn to global bifurcation of solutions of (3.3). By using Theorem 3.3 and Remark 3.4 in Shi and Wang (2009), we conclude that there exists a connected component \(\Lambda ^+\) of \(\Lambda \) containing \(\Gamma _3\), and the closure of \(\Lambda ^+\) includes the bifurcation point \((\delta _*,M_b,0,M_b)\). Moreover, \(\Lambda ^+\) satisfies one of the following three alternatives:

  1. (a)

    it is not compact in \({\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2\);

  2. (b)

    it contains another bifurcation point \(({\hat{\delta }},M_b,0,M_b)\) with \({\hat{\delta }}\ne \delta _*\);

  3. (c)

    it contains a point \((\delta ,M_b+N_3^{**},B_3^{**}(x), M_b+M_3^{**}(x))\) with \(0\ne (N_3^{**},B_3^{**}(x),M_3^{**}(x))\in {\mathcal {W}}\), where \({\mathcal {W}}\) is a closed complement of \(\ker H={{\,\mathrm{span}\,}}({\bar{\zeta }},{\bar{\varphi }},{\bar{\psi }})\) in \({\mathbb {R}}\times X_1\times X_2\).

If the alternative (c) holds, then

$$\begin{aligned} \int _0^{x_h}B_3^{**}(x){\bar{\varphi }}(x)\mathrm{d}x=0. \end{aligned}$$
(A.19)

But it follows from Lemma A.2 that \(B_3^{**}(x)>0\) on \([0,x_h]\), and \({\bar{\varphi }}(x)>0\) as it is a positive eigenfunction. This is a contradiction to (A.19), which means that (c) cannot happen. Suppose that the alternative (b) occurs, then (A.15) has a nonzero solution \(({{\hat{\zeta }}},{{\hat{\varphi }}},{{\hat{\psi }}})\) with \(\delta _*\) replaced by \({\hat{\delta }}\), \({\hat{\varphi \in }} X_1\), and \({{\hat{\varphi }}}>0\). But the eigenvalue problem (3.4) has only one eigenvalue with positive eigenfunction, hence (b) cannot happen either.

Therefore, the alternative (a) must happen, and \(\Lambda ^+\) is not compact in \({\mathbb {R}}^+\times {\mathbb {R}}\times X_1\times X_2\). From Lemma A.2, \((N_3^*,B_3^*(x),M_3^*(x))\) is bounded on \([0,x_h]\) for \(\delta \in [\epsilon ,\delta _{**})\) with any \(\epsilon >0\), and (3.3) has no nonnegative solution when \(\delta >\delta _{**}\). Hence the projection of \(\Lambda ^+\) onto \(\delta \)-axis is contained in \((0,\delta _{**})\), but contains \([\epsilon ,\delta _*)\) for any \(\epsilon >0\), so the projection of \(\Lambda ^+\) onto \(\delta \)-axis contains \((0,\delta _*)\). This also implies that there is at least one positive solution of (3.3) on \(\Lambda ^+\) for any \(\delta \in (0,\delta _*)\). This completes the proof of parts (i) and (ii).

\(\square \)

To prove the existence of coexistence steady state \(E_4\), we first prove the following elementary result.

Lemma A.3

Suppose that \(0<a<x_e(\delta _0^*-\delta _*)\) holds, then

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\varphi ''(x)-v\varphi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _a^*\right) \varphi =0,&{}~0<x<x_h,\\ D_b\varphi '(0)-v\varphi (0)=-a,~D_b\varphi '(x_h)-v\varphi (x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$
(A.20)

has a unique positive solution \({\hat{\varphi }}(x)\).

Proof

By using the transform \(\phi (x)=\varphi (x)e^{-(v/D_b)x}\), we get

$$\begin{aligned} {\left\{ \begin{array}{ll} D_b\phi ''(x)+v\phi '(x)+\left( rf(M_b)g(I(x,0,0))-\delta _a^*\right) \phi =0,&{}~0<x<x_h,\\ -D_b\phi '(0)-a=0,~\phi '(x_h)=0.&{} \end{array}\right. }\nonumber \\ \end{aligned}$$
(A.21)

It is clear that 0 is a lower solution of (A.21). Let \({\tilde{\varphi }}(x)={\tilde{\phi }}(x)e^{(v/D_b)x}\) be the positive eigenfunction of (3.4) corresponding to \(\lambda _1(D_b,v,rf(M_b)g(I(x,0,0))=\delta _*\) with \(\max _{x\in [0,x_h]}{\tilde{\phi }}(x)=1\). Define \({\bar{\phi }}(x)=K [{\tilde{\phi }}(x)+\epsilon (x-x_h)^2]\), where K and \(\epsilon \) are positive constants to be specified. Then

$$\begin{aligned} \begin{aligned}&D_b{\bar{\phi }}''(x)+v{\bar{\phi }}'(x)+\left( rf(M_b)g(I(x,0,0)) -\delta _a^*\right) {\bar{\phi }}\\ =&(\delta _*-\delta _a^*)K{\tilde{\phi }}+\epsilon K \left[ 2D_b+2v(x-x_h)+\left( rf(M_b)g(I(x,0,0)) -\delta _a^*\right) (x-x_h)^2\right] \\\le&0, \;\;\;\; 0<x<x_h, \end{aligned} \end{aligned}$$

if \(\epsilon >0\) is chosen sufficiently small since \(\delta _*-\delta _a^*<0\) and \({\tilde{\phi }}>0\). We can further choose \(K>0\) so that \(-D_b{\bar{\phi }}'(0)-a=2D_bK\epsilon x_h-a=0\) and \({\bar{\phi }}'(x_h)=0.\)

Hence \({\bar{\phi }}\) is an upper solution of (A.21). From Theorem 3.2.1 in Pao (1992), there is a solution \({\hat{\phi }}\) of (A.21) satisfying \(0\le {\hat{\phi }}\le {\bar{\phi }}\). It follows from the strong maximum principle that \({\hat{\phi }}>0\). Note that (A.21) is a linear ODE with nonhomogeneous boundary conditions. Hence \({\hat{\phi }}\) is unique. This implies that there exists a unique positive solution \({\hat{\varphi }}={\hat{\phi }}e^{(v/D_b)x}\) of (A.20). \(\square \)

To obtain the global existence of \(E_4\), we first establish the following a priori estimates for nonnegative solutions \((A_4,N_4,B_4(x),M_4(x))\) of (3.14).

Lemma A.4

Assume that \((A_4,N_4,B_4,M_4)\in ({\mathbb {R}}^+)^2\times Y^2\) is a nonnegative solution of (3.14) with \(A_4,N_4>0\) and \(B_4,M_4\not \equiv 0\). Then

  1. (i)

    \(0<\delta <\delta _a^{**}\), where \(\delta _a^{**}\) is defined in (3.5);

  2. (ii)

    \(0<A_4<K\), where K satisfies \(\displaystyle \frac{r}{x_e}\int _{-x_e}^0g(I(x,K))\mathrm{d}x=\delta +\displaystyle \frac{a}{x_e}\);

  3. (iii)

    For any \(\epsilon >0\), there exists a positive constant \(C(\epsilon )>0\) such that

    $$\begin{aligned}&0<B_4(x)\le C(\epsilon ),\\&0<N_4,M_4(x)<M_b+x_h^2\theta (r+p\delta _a^{**})C(\epsilon )/D_m+ x_h\theta (x_e(1-p)\delta _a^{**}+a)K/D_m \end{aligned}$$

    on \([0,x_h]\) for \(\delta \in [\epsilon ,\delta _a^{**})\).

Proof

It follows from the first equation of (3.14) that (i) and (ii) hold. By applying similar arguments to those in Lemma A.2, we conclude that \(M_4>0\) and for any \(\epsilon >0\), there exists a positive constant \(C(\epsilon )>0\) such that \(0<B_4\le C(\epsilon )\). From the second and fourth equations of (3.14), we have

$$\begin{aligned} b(M_4(0)-N_4)=\theta ( x_e(1-p)\delta +a)A_4 \end{aligned}$$

and

$$\begin{aligned} |D_mM'_4(x)|&=\left| D_m\displaystyle \int _0^{x}M''_4(s)ds+D_mM_4'(0)\right| \\&\le \theta \left| \displaystyle \int _0^x\left( rf(M_4)g(I(x,A_4,B_4))-p\delta \right) B_4ds\right| + b|M_4(0)-N_4|\\&\le x_h\theta (r+p\delta _a^{**})C(\epsilon )+\theta (x_e(1-p)\delta _a^{**}+a)K. \end{aligned}$$

By (A.13), we conclude that (iii) holds. \(\square \)

Proof of Theorem 3.7

Define a nonlinear mapping \(G:{\mathbb {R}}^+\times {\mathbb {R}}^2\times X_3\times X_2\rightarrow {\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3\) by

$$\begin{aligned} G(\delta ,A,N,B,M)=\begin{pmatrix} \left( rf(N)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A))\mathrm{d}x-\delta -\displaystyle \frac{a}{x_e}\right) A\\ \displaystyle \frac{b}{x_e}(M(0)-N)+\theta \left( p \delta -rf(N)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A))\mathrm{d}x\right) A\\ D_bB''-vB'+rBf(M)g(I(x,A,B))-\delta B\\ D_mM''+\theta \left( p\delta - rf(M)g(I(x,A,B))\right) B\\ D_bB'(0)-vB(0)+aA\\ D_mM'(0)-b(M(0)-N)\\ M(x_h)-M_b \end{pmatrix},\nonumber \\ \end{aligned}$$
(A.22)

where \(X_3:=\{u\in C^2([0,x_h]):D_bu'(x_h)-vu(x_h)=0\}\), and \(X_2,Y\) are defined in (3.13). It follows that \(G(\delta ,0,M_b,0,M_b)=0\). Let \(L:=G_{(A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)\) be the Frechét derivative of G with respect to (ANBM). For any \((\xi ,\zeta ,\varphi ,\psi )\in {\mathbb {R}}^2\times X_3\times X_2\), we have

$$\begin{aligned} \begin{aligned} L[\xi ,\zeta ,\varphi ,\psi ] =\begin{pmatrix}0\\ -\theta \left( (1-p)\delta _a^*+\displaystyle \frac{a}{x_e}\right) \xi -\displaystyle \frac{b}{x_e}\zeta +\displaystyle \frac{b}{x_e}\psi (0)\\ D_b\varphi ''(x)-v\varphi '(x)+(rf(M_b)g(I(x,0,0))-\delta _a^*)\varphi \\ D_m\psi ''(x)+\theta (p\delta _a^*-rf(M_b)g(I(x,0,0)))\varphi \\ D_b\varphi '(0)-v\varphi (0)+a\xi \\ D_m\psi '(0)-b(\psi (0)-\zeta )\\ \psi (x_h) \end{pmatrix}. \end{aligned} \end{aligned}$$

If \((\xi _1,\zeta _1,\varphi _1,\psi _1)\in \ker L\), then

$$\begin{aligned} \begin{aligned}&-\theta \left( (1-p)\delta _a^*+\displaystyle \frac{a}{x_e}\right) \xi _1-\displaystyle \frac{b}{x_e}\zeta _1 +\displaystyle \frac{b}{x_e}\psi _1(0)=0,\\&D_b\varphi ''_1(x)-v\varphi '_1(x)+(rf(M_b)g(I(x,0,0))-\delta _a^*)\varphi _1=0,~0<x<x_h,\\&D_m\psi ''_1(x)+\theta (p\delta _a^*-rf(M_b)g(I(x,0,0)))\varphi _1=0,~0<x<x_h,\\&D_b\varphi '_1(0)-v\varphi _1(0)+a\xi _1=0,~D_m\psi '_1(0)-b(\psi _1(0)-\zeta _1)=0,~\psi _1(x_h)=0. \end{aligned} \end{aligned}$$
(A.23)

Let \(\xi _1=1\). Then, from Lemma A.3, \(\varphi _1={\hat{\varphi }}>0\) can be uniquely solved, so are \(\zeta _1={\hat{\zeta }}\) and \(\psi _1={\hat{\psi }}\). Similar to those in the proof of (i) in Theorem 3.5, we have \({\hat{\psi }}\le 0\) and \({\hat{\zeta }}<0\). If \(\xi _1=0\), we can conclude that \(\varphi _1\equiv 0\) as \(\delta _*<\delta _a^*\) and consequently \(\zeta _1=\psi _1=0\). Hence, \(\dim \ker L=1\) and \(\ker L={{\,\mathrm{span}\,}}\{(1,{\hat{\zeta }},{\hat{\varphi }},{\hat{\psi }})\}\). It is also easy to observe that \({{\,\mathrm{codim}\,}}{{\,\mathrm{range}\,}}L=1\) as

$$\begin{aligned} {{\,\mathrm{range}\,}}L=&\left\{ (\rho _1,\rho _2,\rho _3,\rho _4,\rho _5,\rho _6,\rho _7)\in {\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3:\rho _1=0\right\} , \end{aligned}$$

and

$$\begin{aligned} G_{\delta (A,N,B,M)}(\delta _a^*,0,M_b,0,M_b) (1,{\hat{\zeta }},{\hat{\varphi }},{\hat{\psi }}) =(-1,\theta p,-{\hat{\varphi }},\theta p{\hat{\varphi }},0,0,0)\not \in {{\,\mathrm{range}\,}}L. \end{aligned}$$

By applying Theorem 1.7 in Crandall and Rabinowitz (1971), there exists a positive constant \(\varepsilon _4>0\) such that all solutions of (3.14) near \((\delta _*,M_b,0,M_b)\) lie on a smooth curve

$$\begin{aligned} \Gamma _4=\{(\delta _4(s),A_4^*(s),N_4^*(s),B_4^*(s,x),M_4^*(s,x)):0<s<\varepsilon _4\} \end{aligned}$$

with the form

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}A_4^*(s)=s+o(s^2),\;\; N_4^*(s)=M_b+s{\hat{\zeta }}+o(s^2),\\ &{}B_4^*(s,x)=s{\hat{\varphi }}(x)+o(s^2),\;\; M_4^*(s,x)=M_b+s{\hat{\psi }}(x)+o(s^2). \end{array}\right. } \end{aligned}$$
(A.24)

Again the direction of bifurcation can be calculated by

$$\begin{aligned} \begin{aligned} \delta _4'(0)&=-\displaystyle \frac{\left\langle {\hat{l}},G_{(A,N,B,M)(A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)[1,{\hat{\zeta }}, {\hat{\varphi }},{\hat{\psi }}]^2\right\rangle }{2\left\langle {\hat{l}}, G_{\delta (A,N,B,M)}(\delta _a^*,0,M_b,0,M_b)[1,{\hat{\zeta }}, {\hat{\varphi }},{\hat{\psi }}]\right\rangle }\\&=-\displaystyle \frac{r\gamma _2\mu (0)M_b}{M_b+\gamma _1}+\displaystyle \frac{r\gamma _1{\hat{\zeta }}}{(M_b+\gamma _1)^2}\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,0))\mathrm{d}x<0 \end{aligned} \end{aligned}$$

where \({\hat{l}}\) is a linear functional on \({\mathbb {R}}^2\times Y^2\times {\mathbb {R}}^3\) defined as \(\langle {\hat{l}},(\rho _1,\rho _2,\rho _3,\rho _4,\rho _5,\rho _6,\rho _7)\rangle = \rho _1\) and \(\mu (0)\) can been found in (A.8). This shows that the bifurcation at \((\delta _a^*,0,M_b,0,M_b)\) is backward. This proves part (iii). The proof of part (i) and (ii) is similar to the one for Theorem 3.5, so we omit it here. \(\square \)

Appendix B

We briefly describe the numerical algorithm used in the paper. Divide the interval \([0,x_h]\) to n equal size subintervals \([x^i,x^{i+1}]\), \(i=1,2,\ldots ,n,n+1\) with \(x^i=x^{i-1}+\Delta x, x^1=0\) and \(\Delta x=x_h/n\) (grid size). We also denote the step size (in the time direction) by \(\Delta t\). Let

$$\begin{aligned} u_j^i=u(x^i,t_j),~t_j=t_{j-1}+\Delta t,~x^i=x^{i-1}+\Delta x,~1\le j\le m,~1\le i\le n+1. \end{aligned}$$

We also define \(u^0_j\) to be the value at \(x^0=-\Delta x\) and \(u^{n+2}_j\) to be the value at \(x^{n+2}=x_h+\Delta x\). We use backward differences and central differences as

$$\begin{aligned} u_t(x^i,t_j)=\displaystyle \frac{u_{j}^i-u_{j-1}^i}{\Delta t},~u_x(x^i,t_j)=\displaystyle \frac{u_{j}^i-u_{j}^{i-1}}{\Delta x},~ u_{xx}(x^i,t_j)=\displaystyle \frac{u_{j}^{i+1}-2u_{j}^i+u_j^{i-1}}{(\Delta x)^2}. \end{aligned}$$

Then, model (2.1) becomes

$$\begin{aligned} \begin{aligned}&\displaystyle \frac{A_j-A_{j-1}}{\Delta t}=rA_jf(N_j)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x -\displaystyle \frac{a}{x_e}A_j-\delta A_j,\\&\displaystyle \frac{N_j-N_{j-1}}{\Delta t}=\displaystyle \frac{b}{x_e}(M^1_j-N_j)+\theta p\delta A_j-\theta rA_jf(N_j)\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x,\\&\displaystyle \frac{B_{j}^i-B_{j-1}^i}{\Delta t}=D_b\displaystyle \frac{B_{j}^{i+1}-2B_{j}^i+B_j^{i-1}}{(\Delta x)^2}-v\displaystyle \frac{B_{j}^i-B_{j}^{i-1}}{\Delta x}\\&\qquad \qquad \qquad \quad +rB_j^if(M_j^i)g(I(x^i,A_j,B_j^i))-\delta B_j^i,\\&\displaystyle \frac{M_{j}^i-M_{j-1}^i}{\Delta t}=D_m\displaystyle \frac{M_{j}^{i+1}-2M_{j}^i+M_j^{i-1}}{(\Delta x)^2} +\theta p\delta B_j^i-\theta rB_j^if(M_j^i)g(I(x^i,A_j,B_j^i)),\\&D_b\displaystyle \frac{B_{j}^1-B_{j}^{0}}{\Delta x}-vB_j^1=-aA_j,~D_b\displaystyle \frac{B_{j}^{n+2}-B_{j}^{n+1}}{\Delta x}-vB_j^{n+1}=0,\\&D_m\displaystyle \frac{M_{j}^1-M_{j}^{0}}{\Delta x}=b(M_{j}^1-N_j),~D_m\displaystyle \frac{M_{j}^{n+2}-M_{j}^{n+1}}{\Delta x}=M_b,\\&\displaystyle \frac{1}{x_e}\displaystyle \int _{-x_e}^0g(I(x,A_j))\mathrm{d}x=\frac{1}{x_e(K_{bg}+lA_j)}\ln \frac{I_{in}+\gamma _2}{I_{in}\exp (-x_e(K_{bg}+lA_j))+\gamma _2}, \\&I(x^i,A_j,B_j^i)=I_{in}\exp \left( -K_{bg}(x^i+x_e)-x_elA_j-l\sum _{k=0}^i B_j^k\right) , \end{aligned} \end{aligned}$$

for \(1\le j\le m,~1\le i\le n+1\). The simulations in the paper are based on the above numerical algorithm and are implemented in MATLAB.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, J., Kong, J.D., Shi, J. et al. Phytoplankton Competition for Nutrients and Light in a Stratified Lake: A Mathematical Model Connecting Epilimnion and Hypolimnion. J Nonlinear Sci 31, 35 (2021). https://doi.org/10.1007/s00332-021-09693-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1007/s00332-021-09693-6

Keywords

Mathematics Subject Classification

Navigation