Hindawi Publishing Corporation Abstract and Applied Analysis Volume 2014, Article ID 278457, 11 pages http://dx.doi.org/10.1155/2014/278457 Research Article Nonlinear Analysis in a Nutrient-Algae-Zooplankton System with Sinking of Algae ChuanjunDai1,2,3andMinZhao2,3 1KeyLaboratoryofSaline-AlkaliVegetationEcologyRestorationinOilField,MinistryofEducation, AlkaliSoilNaturalEnvironmentalScienceCenter,NortheastForestryUniversity,Harbin150040,China 2ZhejiangProvincialKeyLaboratoryforWaterEnvironmentandMarineBiologicalResourcesProtection, WenzhouUniversity,Wenzhou,Zhejiang325035,China 3SchoolofLifeandEnvironmentalScience,WenzhouUniversity,Wenzhou,Zhejiang325035,China CorrespondenceshouldbeaddressedtoMinZhao;[email protected] Received2April2014;Revised12June2014;Accepted12June2014;Published7July2014 AcademicEditor:MarianoTorrisi Copyright©2014C.DaiandM.Zhao.ThisisanopenaccessarticledistributedundertheCreativeCommonsAttributionLicense, whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited. Areaction-diffusion-advectionmodelisproposedfortheZeyaReservoirtostudyinteractionsbetweenalgaeandzooplankton, includingthediffusivespreadofalgaeandzooplanktonandthesinkingofalgae.Themodelisinvestigatedbothwithandwithout sinking. Conditions of Hopf and Turing bifurcation in the spatial domain are obtained, and conditions for differential-flow instabilitythatgivesrisetotheformationofspatialpatternsarederived.Usingnumericalsimulation,theauthorsexaminethe impactsonalgaeofdifferentnutrientconcentrations,differentsinkingrates,andvariousdiffusivespreadingpatterns.Finally,the modelswithandwithoutsinkingarecompared,revealingthatthesinkingofalgaeplaysanimportantroleintheoscillationsof algaeandzooplankton.AlltheseresultsmayhelptoachieveabetterunderstandingoftheimpactofalgaeintheZeyaReservoir. 1.Introduction Anumberofempiricalandtheoreticalstudiesonpattern formationinreaction-diffusionsystems[7–12]havebeencar- Withtheeconomicdevelopmentofhumansociety,watersin riedoutsincethepioneeringworkofTuring[13].Nowadays, lakesandreservoirsareexperiencingmoreandmoreserious manyapplicationsuseTuring’sscenariotoattempttoexplain eutrophication,whichcancausesustainedalgalgrowth.With patterns in fish skin, vegetation, plankton, and so on [14– a high degree of eutrophication, the rapid growth of algae 16]. However, the reaction-diffusion system cannot explain may form algal blooms, bringing about ecological failure well certain phenomena in vegetation, plankton, and other andevencausingharmtohumans.Forexample,duetothe biologicalsystemsbecausethesesystemsincludeflow.There- eutrophication, algal blooms frequently appear in the Zeya fore, it is necessary and interesting to investigate reaction- Reservoir in Wenzhou, which is located in a subtropical diffusion systems involving flow, which are called reaction- region; this can cause deterioration in water quality, which diffusion-advectionsystems.Systemsexhibitingflowasapart candeprivethedrinkingwaterofmillionsofpeople. ofpatternformation,suchasphytoplankton-nutrientsystems Becauseofthelocalandglobalimpactsofalgaebloomson with sinking, have been recently studied [17]. Therefore, waterquality,thestudiesofdynamicmodelofplanktonare it is practical and significant to study reaction-diffusion- ofcurrentinterest,andnumerousstudieshaveaddressedthe advectionsystems. dynamics of phytoplankton, zooplankton, and algal growth Theorganizationofthispaperisasfollows.InSection2, [1–3]. In a number of these studies, the researchers have wepresentthenutrient-algae-zooplanktonmodel,explaining shown great interest in identifying patterns (or spatiotem- the biological model. In Section3, the model is analyzed, poralcomplexity)ofphytoplanktonandzooplanktongrowth includingtheonewithdiffusionandtheonewithdiffusion [4–6]. and sinking, and the critical conditions of unstability are 2 AbstractandAppliedAnalysis obtained. In addition, a series of simulations are given in Thenutritionalgrowthfunction(𝜇𝑁/(𝐻𝑁+𝑁))ofalgae, Section4. Using simulation, we investigate the effect of as well as the dependence of the zooplankton grazing rate critical factor on the model. Finally, the paper ends with a onalgaedensity,isinagreementwiththeMonodtype[29]. briefdiscussionandconclusionsinSection5. Hence, in the absence of zooplankton, algal growth will saturateat𝐴0 = 𝜇𝑁/[𝑐(𝐻𝑁 +𝑁)],where𝑐 = 𝜇𝑁/[(𝐻𝑁 + 𝑁)𝐾]isthecompetitioncoefficientofthealgaepopulation. 2.Model Growth constraints imposed by different nutrients are not consideredseparately,butanoverallcarryingcapacity,which The reservoir ecology is a relatively ideal ecological envi- depends on the total nutrient level, is assumed [30]. Algae ronment which includes a number of species, mainly phy- predation by zooplankton follows a Michaelis-Menten (or toplankton (with algae composing the majority of the phy- Michaelis-Menten-Holling) type model; in fact, in recent toplankton),zooplankton,andfish.Inareservoir,twofood years,manyexperimentsandobservationshaveshownthat chains typically exist: (I) the grazing food chain: nutrients- the functional response should be described by a ratio- phytoplankton-zooplankton-fish and (II) the detritus food dependentfunction[31,32].Inaddition,manyworksinthe chain: dead organisms-bacteria or benthonic organisms- literature have examined the so-called Michaelis-Menten- zooplankton-fish [18]. These two food chains constitute the typeratio-dependentpredator-preymodel[33–38]. food web, which maintains the ecological balance of the reservoir.Inthisfoodweb,oncethetrophiclevelisabnormal, thebalancemaybedisturbed,forexample,byanalgalbloom. 3.AnalysisoftheModel In recent years, more and more researchers have been interested in studying the nutrients-phytoplankton-zoop- Toperformtheanalysis,model(1)isrewrittenasfollows: lanktonsystembecauseoftheseriousimpactofalgalblooms ontheecologicalenvironmentandonhumanhealth.Anum- 𝜕𝐴 𝜕𝐴 =𝑓(𝐴,𝑍)−V +𝐷 ∇2𝐴, ber of researchers have simulated the relations among phy- 𝜕𝑡 𝜕𝑧 1 toplankton, zooplankton, and fish using reaction-diffusion (2) 𝜕𝑍 models [19–22]. Recently, some studies have included the =𝑔(𝐴,𝑍)+𝐷 ∇2𝑍, sinkingofphytoplanktonintheirmodels[23–25].Sedimen- 𝜕𝑡 2 tationlossesofphytoplanktonseemtobeanimportantfac- tor that may limit phytoplankton biomass at low mixing where𝑓(𝐴,𝑍)=(𝜇𝑁/(𝐻𝑁+𝑁))𝐴−(𝜇𝑁/(𝐻𝑁+𝑁)𝐾)𝐴2− depths [26–28]. Because most phytoplankton taxa have a 𝑎𝐴𝑍/(𝑍+𝑎ℎ𝐴), 𝑔(𝐴,𝑍)=𝜀𝑎𝐴𝑍/(𝑍+𝑎ℎ𝐴)−𝑚𝑍. higherspecificmassthanwater,allnonmotileandnegatively buoyant taxa tend to sink out of the epilimnion with time [23]. The specific sedimentation loss rate of a given taxon 3.1. Model without Diffusion and Sinking. The nonspatial dependsstronglyonitssinkingvelocity,andthereforesinking model has two positive equilibriums under the conditions velocityisanimportantdeterminantofsedimentationlosses; 𝑚ℎ<𝜀and𝜇𝑁𝜀>𝑎(𝐻𝑁+𝑁)(𝜀−𝑚ℎ),whichcorrespondto the actual sedimentation loss rate also depends on mixing spatiallyhomogeneousequilibriumofthefullmodel(1).The depth[23]. equilibriumsofthesystemare Based on the above analysis, consider a vertical water column.Let𝑧indicatethedepthofthewatercolumn.Let𝐴 (i)𝐸1(𝐾,0)(extinctionofthezooplankton); denotethealgaepopulationdensityand𝑍thezooplankton population density. The algae-zooplankton model can be (ii)𝐸∗(𝐴∗,𝑍∗), where 𝐴∗ = 𝐾(𝜇𝑁𝜀 − 𝑎(𝐻𝑁 + 𝑁)(𝜀 − describedbyareaction-diffusionsystemwithadvection: 𝑚ℎ))/𝜇𝑁𝜀,𝑍∗ =(𝑎(𝜀−𝑚ℎ)/𝑚)𝐴∗. 𝜕𝐴 𝜇𝑁 𝜇𝑁 = 𝐴− 𝐴2 The Jacobian matrix, 𝐵, of nonspatial system at the 𝜕𝑡 𝐻 +𝑁 (𝐻 +𝑁)𝐾 𝑁 𝑁 equilibrium𝐸1is 𝑎𝐴𝑍 𝜕𝐴 − −V +𝐷 ∇2𝐴 (1) 𝑍+𝑎ℎ𝐴 𝜕𝑧 1 𝜇𝑁 −1 − 𝑏 𝑏 𝜕𝜕𝑍𝑡 = 𝑍𝜀+𝑎𝐴𝑎𝑍ℎ𝐴 −𝑚𝑍+𝐷2∇2𝑍, 𝐵=(𝑏1211 𝑏1222)=( 𝐻𝑁0+𝑁 ℎ𝜀 −ℎ𝑚). (3) where𝑁isthenutrientlevelofthesystem,𝜇isthespecific growthrate,𝑎isthegrazingrateofzooplanktononalgae,𝜀 From 𝐵, when the condition, 𝑚ℎ > 𝜀, holds, the isthepreyassimilationefficiencyofzooplankton,𝐻𝑁isthe equilibrium𝐸1 islocallystable,whoseindexis+1,andthen nutrientavailabilityconstraint,𝐾isthecarryingcapacity,ℎis there is no positive equilibrium in the nonspatial model. thehandlingtime,𝑚isthemortalityandrespirationrateof When the condition, 𝑚ℎ < 𝜀, holds, the equilibrium 𝐸1 is zooplankton,Visthealgaesinkingvelocity,and𝐷1, 𝐷2 are saddle and then there is a positive equilibrium, 𝐸∗, in the therespectivediffusioncoefficients,∇2 = 𝜕/𝜕𝑥2 +𝜕/𝜕𝑧2.In nonspatialmodelwhenthecondition,𝜇𝑁𝜀>𝑎(𝐻𝑁+𝑁)(𝜀− model(1),allparametersarepositiveconstants. 𝑚ℎ),holds. AbstractandAppliedAnalysis 3 0.05 0.36 Turing I 0.34 0 III 0.32 1 N 𝜆) −0.05 Hopf e ( 0.30 R 2 II 3 −0.10 0.28 IV 4 0.26 −0.15 0.5 1 1.5 2 0 0.3 0.6 0.9 1.2 D2 k (a) (b) Figure1:(a)Bifurcationcurvesofmodel(2)showingtheTuringspace,whichismarkedasIII.(b)Dispersionrelationshipsshowingan unstableHopfmodeandaTuringmode,𝐷2=0.5,with(1)𝑁=0.26, (2)𝑁=0.30,(3)𝑁=𝑁𝑇≈0.33and(4)𝑁=0.36. TheJacobianmatrix,𝐽,ofnonspatialsystemattheequi- where librium𝐸∗is 𝐷 0 𝐷=( 01 𝐷 ). (8) 𝜕 𝑓 𝜕 𝑓 𝑎 𝑎 2 𝐽=( 𝐴 𝑍 ) =( 11 12) 𝜕 𝑔 𝜕 𝑔 𝑎 𝑎 𝐴 𝑍 (𝐴∗,𝑍∗) 21 22 Note that (7) can be solved using the characteristic polynomialoftheoriginalproblem: 𝜇𝑁 𝑎𝑚2ℎ2 −ℎ𝑚2 𝑎− 𝐻 +𝑁 − 𝜀2 𝜀2 (4) 𝜆2−tr𝑘𝜆+Δ𝑘 =0, (9) =( 𝑁 ). 𝑎(𝑚ℎ−𝜀)2 𝑚(𝑚ℎ−𝜀) where 𝜀 𝜀 tr𝑘 =𝑎11+𝑎22−(𝐷1+𝐷2)𝑘2 =tr0−(𝐷1+𝐷2)𝑘2, By 𝐽, the index of 𝐸∗ is +1, and the equilibrium 𝐸∗ is Δ =𝑎 𝑎 −𝑎 𝑎 −(𝑎 𝐷 +𝑎 𝐷 )𝑘2+𝐷 𝐷 𝑘4 (10) locallystablewhenthecondition,𝜇 > max(𝑎(𝐻𝑁+𝑁)(𝜀− 𝑘 11 22 12 21 11 2 22 1 1 2 𝑚ℎ)/𝑁𝜀, (𝑎(𝑚ℎ+𝜀)−𝑚𝜀)(𝐻𝑁+𝑁)(𝜀−𝑚ℎ)/𝑁𝜀2),holds. =Δ0−(𝑎11𝐷2+𝑎22𝐷1)𝑘2+𝐷1𝐷2𝑘4. Therootsof (9)canbeexpressedinthefollowingform: 3.2. Model with Diffusion. Now let us investigate the most ainntaelryessitsi,nmgoeqduelil(ib2)riiusmrepworiinttte𝐸n∗a.sTfoolploewrfos:rmalinearstability 𝜆 (𝑘)= tr𝑘±√tr2𝑘−4Δ𝑘. (11) 1,2 2 𝜕𝐴 =𝑓(𝐴,𝑍)+𝐷 ∇2𝐴, At the bifurcation point, two equilibrial states of the 𝜕𝑡 1 model intersect and exchange their stability. The space- (5) 𝜕𝑍 independentHopfbifurcationbreaksthetemporalsymmetry =𝑔(𝐴,𝑍)+𝐷 ∇2𝑍. 𝜕𝑡 2 of the system, giving rise to oscillations that are uniform in space and periodic in time. In addition, the Turing bifurcation breaks the spatial symmetry, forming patterns Model (5) can be linearizedaroundthespatiallyhomo- geneousfixedpoint𝐸∗forsmallspace-andtime-dependent thatarestationaryintimeandoscillatoryinspace. TheHopfbifurcationoccurswhen fluctuationsandexpandedinFourierspace: Im(𝜆(𝑘)) ≠0, Re(𝜆(𝑘))=0, at 𝑘=0. (12) 𝐴(𝑥⃗,𝑡)∼𝐴∗𝑒𝜆𝑡𝑒𝑖𝑘⃗𝑥⃗, 𝑍(𝑥⃗,𝑡)∼𝑍∗𝑒𝜆𝑡𝑒𝑖𝑘⃗𝑥⃗. (6) Thenthecriticalvalueofthetransition,thecriticalparameter 𝑁onHopfbifurcation,canbeobtainedas Thefollowingcharacteristicequation,then,canbeobtained: 𝐻 (𝑎(𝜀+𝑚ℎ)−𝑚𝜀)(𝜀−𝑚ℎ) 󵄨󵄨󵄨󵄨󵄨𝐽−𝜆𝐼−𝑘2𝐷󵄨󵄨󵄨󵄨󵄨=0, (7) 𝑁𝐻 = 𝜇𝜀2𝑁−(𝑎(𝜀+𝑚ℎ)−𝑚𝜀)(𝜀−𝑚ℎ). (13) 4 AbstractandAppliedAnalysis 6 2 5 0 N = 0.4 4 −2 1 v 3 N = 0.36 00 −4 0 1 𝜆)/ 2 Re ( −6 2 N = 0.32 1 3 −8 0 −10 0 0.05 0.10 0.15 0.20 0.25 k −12 (a) 0 0.05 0.10 0.15 0.20 0.25 0.42 k I 0.40 (a) 2 0.38 0 N 0.36 Spatial patterns −2 1 0.34 −4 0 0 0 0.32 II 1 𝜆)/ −6 2 e ( R 1 2 3 4 5 6 −8 3 (cid:3) −10 (b) Figure3:(a)TypicalneutralcurveV,asdefinedin(27),fordifferent −12 valuesof𝑁.(b)Numericalcalculationofstabilityinthe(V,𝑁)space. Thecurvecanbeobtainedfromthecombinationof(25)and(27). 0 0.05 0.10 0.15 0.20 0.25 k From previous analysis, we can know that the diffusion (b) affectsthestabilityofequilibrium.Next,wewilldiscussthe Figure2:ThedispersionrelationRe(𝜆).Thecriticalmode𝑘 = 𝑘𝑐 effectofsinkingonequilibrium. is obtained for the first unstable point at V = V𝑐 and 𝑁 = 𝑁𝑐, respectively,whileothercurvescorrespondtodifferentvaluesof𝑁 and V. (a) 𝑁 = 0.32, (1)V = 0.2, (2)V = V𝑐, (3)V = 0.1; (b) 3.3.ModelwithDiffusionandSinking. Inthissubsection,the V=0.132,(1)𝑁=0.31, (2)𝑁=𝑁𝑐,(3)𝑁=0.33. modelwithsinkingwillbediscussed.Toanalyzethismodel, followingMurray[15],zero-fluxboundaryconditionswillbe assumedandinitialconditionsspecified.Bysubstituting𝐴= TheTuringbifurcationoccurswhen 𝐴∗+𝐴(𝑟,𝑡)and𝑍=𝑍∗+𝑍(𝑟,𝑡)intomodel(1)andassuming Im(𝜆(𝑘))=0, Re(𝜆(𝑘))=0, at 𝑘=𝑘𝑇≠0. (14) 𝐴 ≪ 1, 𝑍 ≪ 1,alinearapproximationofmodel(1)canbe obtainedasfollows: Thecriticalvalueofthebifurcationparameter𝑁is 𝐻 𝑃 𝜕𝐴 𝜕𝐴 𝑁𝑇 = 𝑄𝑁−𝑃, (15) 𝜕𝑡 =𝑎11𝐴+𝑎12𝑍−V𝜕𝑧 +𝐷1∇2𝐴, (16) w𝐷h1𝐷er2e𝜀𝑃2𝑘=𝑇4,((𝑄𝐷=2𝑎𝜇(𝜀𝜀+(𝑚𝑚(ℎ𝜀)−−𝑚𝜀𝑚ℎ)𝐷+1)𝜀𝑘𝐷𝑇22+𝑘𝑚𝑇2)𝑎,(𝑘𝜀𝑇2−=𝑚√ℎΔ))(0𝜀/𝐷−𝑚1𝐷ℎ2).− 𝜕𝜕𝑍𝑡 =𝑎21𝐴+𝑎22𝑍+𝐷2∇2𝑍. AbstractandAppliedAnalysis 5 0.085 0.3 0.0845 0.25 0.084 0.2 0.0835 0.083 0.15 0.0825 0.1 0.082 0.05 0.0815 (a) (b) 0.2 0.15 0.1 0.05 (c) Figure4:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁=0.25, 𝐷2 =0.5.(a)0iterations, (b)50000iterations,and(c)100000iterations. Theinitialconditionsareassumedasfollows:𝐴|𝑡=0 =𝜑(𝑟) whereΦ(𝑞)andΨ(𝑞)arethetransformsof𝜑(𝑟)and𝜙(𝑟).To and 𝑍|𝑡=0 = 𝜙(𝑟), where the functions 𝜑(𝑟) and 𝜙(𝑟) decay revealthepresenceofinstabilityanddiscloseitsnature,itis rapidly for 𝑟 → ±∞. Following the standard approach, sufficienttoconsideronevariable.Furthercalculationsyield Laplace transformation of the linearized equations over the 𝐴(𝑟,𝑡) twoindependentvariables𝑟and𝑡isperformed.For𝑟,theso- called two-sided version of the transformationis used [39]. 1 𝛼+𝑖∞ =− ∫ 𝑒𝜆𝑡𝑑𝜆 Therelationsfortheforwardandbackwardtransformsare 4𝜋2 𝛼−𝑖∞ ∞ ∞ (20) 𝐴𝜆𝑞 =∫ 𝑒−𝜆𝑡𝑑𝑡∫ 𝐴(𝑟,𝑡)𝑒−𝑞𝑟𝑑𝑞, (17) ×∫𝑖∞ 1 [(𝜆−𝑎 −𝑑 𝑞2)Φ(𝑞) 0 −∞ 𝐷(𝜆,𝑞) 22 2 −𝑖∞ 1 𝛼+𝑖∞ 𝑖∞ 𝐴(𝑟,𝑡)=− ∫ 𝑒𝜆𝑡𝑑𝑠∫ 𝐴 𝑒𝑞𝑟𝑑𝑞, (18) +𝑎 Ψ(𝑞)]𝑒𝑞𝑟𝑑𝑞, 4𝜋2 𝜆𝑞 12 𝛼−𝑖∞ −𝑖∞ wherethedenominatoris where𝜆and𝑞arecomplexvariables.In(17)forthebackward 𝐷(𝜆,𝑞)=(𝜆−𝑎 −𝐷 𝑞2)(𝜆−𝑎 +V𝑞−𝐷 𝑞2)−𝑎 𝑎 . transformation,theintegrationcontourinthe𝑞-planeisthe 22 2 11 1 12 21 imaginaryaxis.Inthe𝜆-plane,thecontourisparalleltothe (21) imaginaryaxisandlocatedtotherightofallsingularitiesof Hence,thedispersionequationis theintegrand. After this transformation, the kinetic equations read as (𝜆−𝑎22−𝐷2𝑞2)(𝜆−𝑎11+V𝑞−𝐷1𝑞2)−𝑎12𝑎21 =0. (22) follows: For the sake of convenience, the one-dimensional case is (𝜆−𝑎 −𝐷 𝑞2+V𝑞)𝐴 −𝑎 𝑍 =Φ(𝑞), consideredfirst.Fromthequadraticequation(22),theroots 11 1 𝜆𝑞 12 𝜆𝑞 (19) 1 (𝜆−𝑎22−𝐷2𝑞2)𝑍𝜆𝑞−𝑎21𝐴𝜆𝑞 =Ψ(𝑞), 𝜆= 2(tr𝑘+V𝑘𝑖±√𝑀+𝑄𝑖) (23) 6 AbstractandAppliedAnalysis 0.1495 0.16 0.149 0.155 0.1485 0.15 0.148 0.145 0.1475 0.14 0.147 0.135 0.1465 0.13 0.146 (a) (b) 0.25 0.2 0.15 0.1 0.05 (c) Figure5:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁=0.31, 𝐷2 =0.5.(a)0iterations, (b)50000iterations,and(c)100000iterations. canbeobtained,where for any values of 𝑘≠0, then pattern formation will occur. Therefore, it is necessary to find the critical value of 𝑘 for 𝑀= 2(𝑎 −𝑎 )(𝐷 −𝐷 )𝑘2 11 22 2 1 Re(𝜆) = 0. In fact, what is needed is the maximum value ofRe(𝜆).Thus,thecriticalconditioncanbeobtainedusing +(𝐷12+𝐷22)𝑘4−2𝐷1𝐷2𝑘2 Re(𝜆)=0,asfollows: (24) +tr20−4Δ0−V2𝑘2, V2 = tr4𝑘−tr2𝑘𝐿, (27) 𝑃2−tr2𝑘2 𝑄=2((𝑎 −𝑎 )𝑘+(𝐷 −𝐷 )𝑘3)V. 𝑘 22 11 2 1 where Straightforwardmanipulationof (23)yields 𝐿= 2(𝑎 −𝑎 )(𝐷 −𝐷 )𝑘2+(𝐷2+𝐷2)𝑘4 Re(𝜆)= 21(tr𝑘±√12(√𝑀2+𝑄2+𝑀)), (25) 11 22 2 1 1 2 −2𝐷1𝐷2𝑘2+tr20−4Δ0, (28) 1 1 Im(𝜆)= (V𝑘±sign(𝑀)(√ (√𝑀2+𝑄2−𝑀))), 𝑃=(𝑎 −𝑎 )𝑘+(𝐷 −𝐷 )𝑘3. 2 2 22 11 2 1 (26) Generically, (27) has a strictly positive minimum V𝑐. whereRe(𝜆)andIm(𝜆)denotetherealandimaginaryparts When (27) holds, the neutral curve V = V(𝑘,𝑁) has a of𝜆,respectively. minimumatanonzerovalueofV=V𝑐 >0. The condition for a spatial mode to be unstable is that From previous analysis, we can know that the sinking model(1)growsintoapattern;thatis,Re(𝜆)>0.IfRe(𝜆)>0 affectsthestabilityofequilibrium. AbstractandAppliedAnalysis 7 0.15 0.25 0.149 0.2 0.148 0.15 0.1 0.147 0.05 0.146 (a) (b) 0.25 0.2 0.15 0.1 0.05 (c) Figure6:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁=0.31,𝐷2 =1.25.(a)0iterations, (b)50000iterations,and(c)100000iterations. 4.NumericalResults biologicalpointofview.Intherestofthispaper,themodel parametersaresetto In the previous, we analyzed the model, including nonspa- 𝜇=0.6, 𝐻 =1, ℎ=2, 𝑎=0.5, tial model, the model with diffusion, and the model with 𝑁 (29) diffusion and sinking. In this section, in order to carry out 𝜀=0.5, 𝑚=0.2, 𝐾=0.5, 𝐷 =0.05. numericalsimulations,itisnecessarytodiscretizethespace 1 andtimevariablesintheproblem.So-called“discretization” Firstly,weconsidertheHopfandTuringbifurcationsof means that the continuous problem is transformed into a model(1)withoutsinkingintheparameterspacespannedby discrete problem; that is, the continuous space becomes a the parameters 𝑁 and 𝐷2, as illustrated in Figure1(a). The discretedomainwith𝑚×𝑛latticesites.Thespacingbetween HopfbifurcationlineandtheTuringbifurcationlineintersect the lattice points is defined by the lattice constant Δℎ. The atacodimension-2bifurcationpoint,calledtheTuring-Hopf following discussion uses an explicit Euler method with bifurcation point. The parametric space can be separated time step Δ𝑡 = 0.01 for time integration and a finite- intofourdomains.DomainI,locatedabovebothbifurcation difference method for the diffusion and advection terms, lines,correspondstosystemswithhomogeneousequilibrium in which the difference approaches the derivatives when which is unconditionally stable. Domain II is the region of Δℎ → 0. The numerical simulation is performed over pureHopfinstability,anddomainIIIisthatofpureTuring 200 × 200 mesh points with spatial resolution Δℎ = 0.5, instability.InDomainIV,bothHopfandTuringinstabilities and the zero-flux boundary condition is used. Because the occur. Figure1(b) shows the dispersion relationships of the spatiotemporaldynamicsofthesystemdependonthechoice unstableHopfmodeandtheTuringmodeastheymovefrom of initial conditions, random spatial distributions of the stabletounstablestates,usingmodel(2). speciesaroundthenontrivialstationarystate𝐸∗werechosen, By (23), the critical wavelength (𝑘𝑐 = 0.089) can be which seems to be a fairly general assumption from the obtained such that Re(𝜆(𝑘 = 𝑘𝑐)) = 0, where 𝑁 and V 8 AbstractandAppliedAnalysis 0.3 0.085 0.25 0.2 0.084 0.15 0.083 0.1 0.082 0.05 0 (a) (b) 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 (c) Figure7:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁 = 0.25, 𝐷2 = 0.5, V = 0.1.(a)0 iterations,(b)50000iterations,and(c)100000iterations. are chosen as control parameters, and their corresponding becausethemodelhasarelativelylargenumberofparame- valuesare𝑁𝑐 = 0.32andV𝑐 = 0.132,asshowninFigure2. ters.However,thenumericalsimulationsshowthatthetype Figure2depictstherangeofvaluesof𝑘forcertainparameter ofdynamicsinthesystemisdeterminedbythevaluesof,𝐷2, values,anditisapparentthatspatialpatternscanoccurwhen andV,whichrepresenttheeffectofnutrients,thediffusionof Re(𝜆) > 0.Atthesametime,itiseasytoseethatinstability zooplankton,andthesinkingofalgae. occurswhenV>V𝑐(Figure2(a)),whileinstabilityalsooccurs Figure4showsthestripelikepatternsthatspontaneously when𝑁<𝑁𝑐(Figure2(b)). form at 𝑁 = 0.25, 𝐷2 = 0.5. From Figure4, it is According to (27), the neutral curve V = V(𝑘,𝑁) has a clearthestripelikespatialpatternsarisefromrandominitial minimum at a nonzero value of V = V𝑐 > 0, as shown in conditions.Aftersometime,thestripelikepatternsformfrom Figure3(a).FromFigure3(a),itiseasytoseethatV𝑐increases the initial state (Figure4(b)) and grow steadily over time with increasing nutrient concentration 𝑁. The reason for untiltheyreachacertainwidth.Finally,thestripelikespatial this is that the algae may pursue nutrients that have been patternsprevailoverthewholedomain(Figure4(c)). depositedonthebottom.Inthisway,nutrientscaninfluence Figure5 shows snapshots of algae spatial patterns with thesinkingofalgae.Moreover,using(27),aparameterspace 𝑁 = 0.31, 𝐷2 = 0.5 at 0, 50000, and 100000 iterations. ((V,𝑁) space) can be obtained in which the stable and From Figure5, it is interesting to note that the spotted unstabledomainsofmodel(1)canbedetermined,asshown spatialpatternsarisefromtherandominitialconditions.The inFigure3(b).Figure3(b)alsoshowsthatV𝑐 increaseswith stripelike pattern forms from the initial state (Figure5(a)) increasing nutrient concentration 𝑁. Moreover, instability after some time has passed. Surprisingly, however, the final will occur for any 𝑁 < 𝑁𝑐 when V𝑐 is fixed, as is shown in pattern becomes a spotted spatial pattern. Compared with Figure3(b)asdomainII(identifiedasthe“spatialpatterns” Figure4, only the amount of nutrients has been changed, domain). butthisisclearlyanessentialdifference.Therefore,nutrient Itisobviousthatdetailednumericalinvestigationofthe concentrationevidentlyplaysanimportantroleinthespatial model in the whole parameter space is virtually impossible distributionofalgae. AbstractandAppliedAnalysis 9 0.12 0.085 0.11 0.1 0.084 0.09 0.083 0.08 0.07 0.082 0.06 0.05 (a) (b) 0.25 0.2 0.15 0.1 0.05 (c) Figure8:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁 = 0.25, 𝐷2 = 0.5, V = 2.(a)0 iterations,(b)50000iterations,and(c)100000iterations. Figure6 shows snapshots of algae spatial patterns with and the spatial distribution of algae. Moreover, the sinking 𝑁 = 0.31, 𝐷2 = 1.25 at 0, 50000, and 100000 iterations. of algae can play an important role in the oscillations of Comparing Figures 5 and 6, the final pattern in Figure6 algae and zooplankton. Such evolutionary developments of is similar to that in Figure5 (they are both spotted), but the reaction-diffusion-advection system may well be taking theirformativeprocessesaredifferent.Moreover,itisobvious placeintherealworld. that the density of spots is sparser in Figure6(c) than in Figure5(c). In fact, the difference between Figures 5 and 6 isthechangeindiffusioncoefficient.Theseresultsshowthat 5.DiscussionandConclusions thediffusionofzooplanktonisakeyfactorintheformative Inthepresentpaper,weusedareaction-diffusion-advection processesofalgaepatterns. model to investigate the interaction between algae and Now including sinking in model (1), that is, V≠0, the zooplankton, where the sinking of algae was considered effects of sinking are simulated, and the results are shown which was described by the advection term. Using linear inFigures7and9.Figure7showssnapshotsofalgaespatial analysis technique, the model was analyzed, and the effects patterns with 𝑁 = 0.25, 𝐷2 = 0.5, V = 0.1 at 0, 50000, ofcriticalfactors,including𝑁, 𝐷2,andV,onthemodelwere and100000iterations.InFigure7,verticalstripesmovefrom discussed.ThentheinstabilityconditionsofHopfandTuring toptobottomandthenexitfromthebottom.Comparedto were obtained. In addition, from the analysis described in Figure4,itisapparentthatsinkingexistsinFigure7.Tosee Section3 and Figures 1 and 9, it was not difficult to find the effect of sinking, V = 0.1 is replaced by V = 2, and thatthetheoreticalresultswereconsistentwiththenumerical the resulting patterns are shown in Figure8. If Figure9 is simulations. comparedwithFigure6,thespottedpatternbecomesvertical The model was relatively simple because it was only stripes.NotethatFigures7,8,and9arequitesimilaratthe an abstraction of real-world phenomena, but it reproduced end. many features of real-world phenomena. In the real world, Basedontheaboveanalysis,itisapparentthatparameters many patterns have been observed, such as banded vegeta- 𝑁, 𝐷2,andVarecriticalfactorsforcontrollingtheamount tion, patches, and spiral waves, which could be regular or 10 AbstractandAppliedAnalysis 0.15 0.1485 0.149 0.148 0.148 0.1475 0.147 0.147 0.146 (a) (b) 0.25 0.2 0.15 0.1 0.05 (c) Figure9:Snapshotsofcontourdiagramsofthetimeevolutionofalgaeatdifferentpointsintimewith𝑁=0.31, 𝐷2 =1.25, V=2.5.(a)0 iterations,(b)1500iterations,and(c)100000iterations. irregular. The patterns may be forced to occur by physical ConflictofInterests factorsorinternalfactors.Whileourexplanationfocusedon apredator-preyinteractionbetweenalgaeandzooplankton, The authors declare that there is no conflict of interests especially, how the sinking of algae affected the interaction regardingthepublicationofthispaper. and patterns, our analytical results showed that the homo- geneous steady state became unstable due to the sinking Acknowledgments of algae. In addition, diffusion also led to instability of the homogeneous steady state. But, the sinking of algae forced This work was supported by the National Natural Science band patterns to occur. That is, when the sinking flux of Foundation of China (Grant no. 31170338), by the Key algae was beyond a certain critical value which could be ProgramofZhejiangProvincialNaturalScienceFoundation obtained from our analyzed results, the spatial distribution of China (Grant no. LZ12C03001), and by the National Key ofpopulationwasnothomogeneousbecausethesinkingof Basic Research Program of China (973 Program, Grant no. algaeexisted. 2012CB426510). Fromthenumericalresults,itwasobviousthatdifferent 𝑁, 𝐷2, or V can lead to different patterns at the same time. Manyfieldstudiesindicatedthattheoscillationofpopulation References density occurs in reality, while Figures4and9 implied that nutrient, diffusion, and the sinking can result in oscillation [1] H.Malchow,“Spatio-temporalpatternformationinnonlinear of population density. It means that these factors may play non-equilibriumplanktondynamics,”ProceedingsoftheRoyal animportantroleinthechangeofalgaepopulationdensity. SocietyB:BiologicalSciences,vol.251,no.1331,pp.103–109,1993. Theseresultsmayprovidebetterunderstandingonthestudy [2] H.Malchow,F.M.Hilker,S.V.Petrovskii,andK.Brauer,“Oscil- ofalgaldynamicsinwaterenvironment. lationsandwavesinavirallyinfectedplanktonsystem—partI: