NPZ-model with seasonal variability of plankton population
dynamics
Citation for published version (APA):
Goorden, S. A., Korzilius, S. P., Thije Boonkkamp, ten, J. H. M., Anthonissen, M. J. H., & Rathish Kumar, B. V. (2011). NPZ-model with seasonal variability of plankton population dynamics. (CASA-report; Vol. 1130). Technische Universiteit Eindhoven.
Document status and date: Published: 01/01/2011
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne Take down policy
If you believe that this document breaches copyright please contact us at: openaccess@tue.nl
providing details and we will investigate your claim.
EINDHOVEN UNIVERSITY OF TECHNOLOGY
Department of Mathematics and Computer Science
CASA-Report 11-30
April 2011
NPZ-model with seasonal variability of plankton population dynamics
by
S.A. Goorden, S. P. Korzilius, J.H.M. ten Thije Boonkkamp,
M.J.H. Anthonissen, B.V. Rathish Kumar
Centre for Analysis, Scientific computing and Applications
Department of Mathematics and Computer Science
Eindhoven University of Technology
P.O. Box 513
5600 MB Eindhoven, The Netherlands
ISSN: 0926-4507
NPZ-model with seasonal variability of plankton
population dynamics
S.A. Goorden1, S.P. Korzilius1, J.H.M. ten Thije Boonkkamp1, M.J.H. Anthonissen1, B.V. Rathish Kumar2
1Department of Mathematics and Computer Science, Eindhoven University of Technology
P.O. Box 513, 5600 MB Eindhoven, The Netherlands
2Department of Mathematics and Statistics, Indian Institute of Technology Kanpur
Kanpur-208016, India
Abstract
In this paper we include seasonal variability in the NPZ-model of a population of phytoplankton. For the space discretization of the governing equations we employ the complete flux scheme. A stability analysis of the fully discrete scheme is presented. A new boundary condition, based on the complete flux scheme, is introduced. Numerical results show good agreement with previous results in literature.
Keywords. NPZ-model, seasonal variability, complete flux scheme, stability analysis, boundary treatment.
1
Introduction
The topic of phytoplankton has been widely studied for the last decades. Phytoplankton are small plants which inhabit the oceans and lakes of the earth. They are too small to be seen by the naked eye, but they cause a green discolouration of the water when present in large amounts. Phytoplankton are important for several reasons. First, they can extract carbon dioxide from the atmosphere and produce oxygen, through the process of photosynthesis. In fact, phytoplankton consumes as much carbon dioxide as all land vegetation does together. Therefore, it has a major influence on climate change [13]. Second, phytoplankton are also important because they are at the bottom of the food chain. Phytoplankton is eaten by marine life, ranging from microscopic organisms to whales. Small fish feed on phytoplankton and they are in turn eaten by larger fish. Some species of phytoplankton can produce biotoxins, which might end up in fish caught for human consumption. Finally, after a bloom, large amounts of dead phytoplankton can sink to the bottom of the ocean or lake. Bacteria that decompose the phytoplankton can extract so much oxygen from the water that it suffocates marine life [10]. All these factors make it important to understand how a population of phytoplankton behaves under different circumstances.
For the concentration of phytoplankton we are only interested in the mixed layer, which stretches from the water surface to a depth of approximately 100 meters. The reason for this is that the water below this layer does not receive enough sunlight for photosynthesis, which phytoplankton need to survive. Phytoplankton leave the mixed layer by sinking slowly. This sinking is descibed by downward advection. Besides sunlight phytoplankton also need nutrients, such as iron, to survive. These nutrients sojourn in the deeper waters of the ocean and enter the mixed layer due to mixing. During summer, when the top
2 MATHEMATICAL MODEL 2
layer is warmer than the water below, there will be less mixing than in winter, when the temperature difference is smaller. Mixing of water is described by diffusion. Another effect of warmer water is that phytoplankton reproduces faster in summer than in winter. Because nutrients are so light, and therefore hardly sink in comparison with phytoplankton, there is no advection term for nutrients. In addition to phytoplankton and nutrients we also include zooplankton in our model, which are microscopic organisms that feed on phytoplankton. They do not need sunlight and do not sink because they have some capability of swimming, i.e., they are not advected downward.
All governing equations are of (advection-)diffusion-reaction type. Previous research focussed on various aspects, such as stability of the governing PDEs [19], inclusion of seasonal variations [17], one-dimensional models [14] or multi-one-dimensional models [13]. In this paper we adopt a one-one-dimensional model, which includes several features of other models reported in literature, and in particular seasonal variability. The numerical scheme we use is the finite volume-complete flux scheme, described in [15, 16]. The complete flux scheme is a flux approximation scheme for the advection-diffusion-reaction equation. The basic idea is to compute the numerical flux from a local BVP for the entire equation, including the source term. As a result, the numerical flux consists of two part, i.e., the homogeneous and inhomogeneous flux, corresponding to the advection-diffusion operator and the source term, respectively. The (constant coefficient) homogeneous flux is used in many applications as numerical flux; see e.g. [11]. The inhomogeneous flux is an extension, and is particulary of importance for dominant advection in combination with a strong source term. It is therefore anticipated that the scheme will work well for turbulence models of phytoplanton, where the (horizontal) advection is dominant over diffusion [7].
In Section 2 we describe our mathematical model, and subsequently in Section 3, we outline the discretisation scheme. A stability analysis of the scheme is presented in Section 4. In Section 5 we introduce the boundary conditions. Numerical results are given in Section 6, and finally, concluding remarks are formulated in Section 7.
2
Mathematical model
In our model we incorporate phytoplankton, zooplankton and nutrients. Phytoplankton consumes nutri-ents and zooplankton feeds on phytoplankton. Phytoplankton needs light for photosynthesis and there-fore it can only live near the surface of the water. Because of gravity, phytoplankton tends to sink slowly, however, near the surface it can be mixed upwards with the water due to wind stress. Plank-ton mostly lives in this mixed layer, which is generally not deeper than 100 meters. We assume that the (hydro)dynamics of the mixed layer does not depend on the horizontal directions, but only on the vertical space coordinate x. The key assumption in this article is that the dynamics of the mixed layer changes with the seasons. In summers the water is warmer, which has two effects. First, phytoplankton reproduces faster than in winter, see e.g. [4], and second, because of higher temperature gradients, there is less mixing of the water, see e.g. [1] and [14].
The governing equations for the concentrations of phytoplankton P [mmol N2m−3], zooplankton Z
[mmol N2m−3] and nutrients N [mmol N2m−3] are given by:
∂P ∂t + w ∂P ∂x = ∂ ∂x K(x, t)∂P ∂x + p(x, t, N ) − `PP − g(P )Z, (2.1a) ∂Z ∂t = ∂ ∂x K(x, t)∂Z ∂x − `ZZ + γg(P )Z, (2.1b) ∂N ∂t = ∂ ∂x K(x, t)∂N ∂x − αp(x, t, N )P + rα `PP + `ZZ, (2.1c)
2 MATHEMATICAL MODEL 3
where w > 0 [m s−1] is the constant sinking velocity of phytoplankton, K(x, t) [m2 s−1] is the eddy diffusivity coefficient, p(x, t, N ) [s−1] is the production rate of phytoplankton and g(P ) [s−1] is the grazing rate, i.e., the rate at which zooplankton eats phytoplankton. Moreover, system (2.1) contains the following coefficients: the mortality rates `P [s−1] and `Z [s−1] of phytoplankton and zooplankton,
respectively, the assimilation coefficient γ [-], which is the fraction of phytoplankton that is converted into new zooplankton, the conversion factor α [-], which is the amount of nutrients that is needed to produce phytoplankton and the regeneration rate r [-], i.e., the fraction of phytoplankton and zooplankton that is converted into nutrients. The advection term for the zooplankton is absent, because it has some capability of staying near the surface. There is also no advection term for nutrients, because they do not have enough mass to sink with a significant velocity. The system (2.1) is referred to as the NPZ-model.
Next, we introduce models for K, p and g. Mixing is generally described by a diffusion term, with eddy diffusivity coefficient K given by
K(x, t) = Ks(x) 1 + µKsin 2πtT , (2.2a) Ks(x) = Kmax 1 + C2 exp − x−x0 C1 2 + C2, (2.2b)
where Ks(x) [m2s−1] is the stationary eddy diffusivity. Parameters in (2.2) are: the diffusivity variability
coefficient µK[-], the period of seasonal variation T , corresponding to one year, the maximum diffusivity
Kmax[m2 s−1] occuring at x = x0, C1[m], which determines the thickness of the mixed layer and C2
[m], which determines the diffusivity at the bottom of the mixed layer. The diffusivity depends on many things, such as wind stress and temperature of the water. We are not modelling any particular geographical location and therefore we have chosen a representative profile for Ks(x), see Figure 1,
which is based on profiles found in literature, see e.g. [6] and [14]. The sinusoidal term in (2.2a) is added to take into account that there is less mixing in summers than in winters. Assuming that we consider the Northern Hemisphere, the water temperature is maximal around August 1 and it is minimal around February 1. This means that, with our definition of K, the initial time t = 0 corresponds to November 1. 0 20 40 60 80 100 0 0.5 1 1.5x 10 −3 x [m] Ks (x) [m 2 s −1 ]
2 MATHEMATICAL MODEL 4
Analogous to the expression for K, the production rate p exhibits a seasonal variation given by
p(x, t, N ) = ps(x, N ) 1 − µpsin 2πtT , (2.3a) ps(x, N ) = pmaxmin Line−νx LH+ Line−νx , N NH+ N , (2.3b)
where ps(x, N ) [s−1] is the stationary production rate. Furthermore, µp[-] is the production variability
coefficient, pmax[s−1] is the maximum production rate, Lin[mol photons m−2s−1] is the incident light
intensity, ν [m−1] is the light absorption coefficient and LH [mol photons m−2 s−1] and NH [mmol
N2 m−3] are the half saturation constants of light and nutrients, respectively. It is assumed that the
production rate of phytoplankton is limited by the amount of light and nutrients available. Both are described according to the Holling type II functional response. For the nutrients this is widely used, see e.g. [1, 14]. For the dependence on light there are several options, and we choose the same Holling type II functional response as in [12]. The light intensity is assumed to decrease exponentially with depth. Some limiting profiles of the stationary production rate are presented in Figure 2. The sinusoidal term in (2.3a) is included to account for the higher production rate in summers.
0 0.02 0.04 0.06 0.08 0.1 0 0.5 1 1.5 2x 10 −5 N [mmol N2 m−3] ps [s −1 ]
(a) unlimited light
0 20 40 60 80 100 0 0.5 1 1.5 2x 10 −5 x [m] ps [s −1 ] (b) unlimited nutrients
Figure 2: Profiles of the stationary production rate ps(x, N ) in case of unlimited light and nutrients,
respectively.
Finally, the grazing rate is defined as
g(P ) = G max 1 − exp − kiv(P − Pmin), 0, (2.4)
where G [s−1] is the phytoplankton predation rate, kiv
mmol N2
−1
m3 is Ivlev’s coefficient and Pmin
[mmol N2m−3] is the treshold concentration of phytoplankton. This term has been taken from [14]. As
one might expect, the grazing rate approaches G when there is sufficient phytoplankton and reduces to 0 when there is very little phytoplankton. A typical profile is given in Figure 3. All parameter values in this section are listed in Table 1.
3 THE COMPLETE FLUX SCHEME 5 0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 x 10−6 P [mmol N2 m−3] g(P) [s −1 ]
Figure 3: Profile of the grazing rate g(P ).
3
The complete flux scheme
In this section we outline the derivation of the complete flux scheme for (2.1); for a detailed account see [15, 16]. We first present the stationary flux approximation and subsequently its extension to time-dependent problems.
Consider the stationary model equation d dx wϕ − K(x)dϕ dx = s, (3.1)
for ϕ = P, Z or N and where s is the corresponding source term from (2.1). Introducing the flux f = wϕ − K(x)dϕ
dx, (3.2)
we can write df /dx = s. We cover the domain with a finite number of disjunct cells (intervals) Ij =
xj−1/2, xj+1/2 and define grid points xj in the centre of each interval where the variable ϕ has to
be approximated. Integrating equation (3.1) over the interval Ij and adopting the midpoint rule for the
integral of the source term, we obtain the discrete equation
Fj+1/2− Fj−1/2= sj∆x, (3.3)
where Fj+1/2 is the numerical flux approximating f at the cell boundary x = xj+1/2 and where sj =
s(xj). We require that Fj+1/2depends on both ϕ and s in the adjacent grid points xj and xj+1, i.e., we
are looking for an expression of the form
Fj+1/2= αj+1/2ϕj− βj+1/2ϕj+1+ ∆x γj+1/2sj+ δj+1/2sj+1, (3.4)
where the coefficients αj+1/2 etc. depend on w and K. A similar expression holds for Fj−1/2.
Re-lation (3.4) is the numerical approximation of the integral representation of f (xj+1/2), which can be
derived from the local BVP on (xj, xj+1) for the entire equation (3.1), including the source term.
Con-sequently, the numerical flux Fj+1/2 is the superposition of the homogeneous flux Fj+1/2h , depending
3 THE COMPLETE FLUX SCHEME 6
of the source term. Introducing the function R(x) = w∆x/K(x) as generalisation of the cell-Reynolds number, we obtain the following expressions
Fj+1/2= Fj+1/2h + Fj+1/2i , (3.5a) Fj+1/2h = − ˜ Rj+1/2 ¯ Rj+1/2 ˜ Kj+1/2 ∆x B ¯Rj+1/2ϕj+1− B − ¯Rj+1/2ϕj , (3.5b) Fj+1/2i = 1 2− W ¯Rj+1/2∆x sj, (3.5c)
where B(z) := z/ ez− 1, W (z) := ez− 1 − z/ z ez− 1, see Figure 4, and where the average ¯
vj+1/2and weighted average ˜vj+1/2(v = R, K) are defined as follows
¯
vj+1/2 = 12 vj+ vj+1, ˜vj+1/2= W − ¯Rj+1/2vj + W ¯Rj+1/2vj+1. (3.6)
Note that the weighted average ˜vj+1/2reduces to the arithmetic average ¯vj+1/2when w = 0 (no advec-tion) and to the upwind value vjwhen K(x) = 0 (no diffusion). We refer to the flux approximation (3.5)
as the complete flux scheme.
Combining (3.3) and (3.5) we obtain the following discretisation scheme
−aW,jϕj−1+ aC,jϕj − aE,jϕj+1= bW,jsj−1+ bC,jsj, (3.7a)
with coefficients aW,jetc. given by
aW,j = ˜ Kj−1/2 ∆x ˜ Rj−1/2 ¯ Rj−1/2 B − ¯Rj−1/2, aE,j = ˜ Kj+1/2 ∆x ˜ Rj+1/2 ¯ Rj+1/2 B ¯Rj+1/2, aC,j = aW,j+1+ aE,j−1, bW,j = ∆x 12 − W ¯Rj−1/2, bC,j = ∆x 12+ W ¯Rj+1/2. (3.7b) The discretisation scheme in (3.7a) gives rise to a linear system of the form
Aϕ = Bs + b, (3.8) −10 −8 −6 −4 −2 0 2 4 6 8 10 2 4 6 8 10 −10 −8 −6 −4 −2 0 2 4 6 8 10 0.2 0.4 0.6 0.8 1
4 STABILITY ANALYSIS 7
with A = tridiag aW,j, aC,j, aE,j and B = tridiag bW,j, bC,j, 0. The matrix A is an M-matrix for
the special case K(x) = Const.
Next, consider the time-dependent extension of (3.1), i.e., ∂ϕ ∂t + ∂ ∂x wϕ − K(x, t)∂ϕ ∂x = s, (3.9)
which can be written as ∂ϕ/∂t + ∂f /∂x = s with f defined in (3.2), however with a time dependent diffusivity coefficient K(x, t). Analogous to the derivation of (3.3), we integrate equation (3.9) over the cell Ij to obtain the semidiscretisation
∆x ˙ϕj+ Fj+1/2− Fj−1/2= sj∆x, (3.10)
where ˙ϕj(t) denotes the approximation of ∂ϕ/∂t(xj, t). The procedure to derive the numerical flux is
completely analogous to the stationary case, except that we include the time derivative in the source term, i.e., we replace s by ˜s = s − ∂ϕ/∂t and determine the numerical flux from a local BVP on (xj, xj+1)
for the quasi-steady equation ∂x wϕ − K(x, t)∂xϕ = ˜s. This way we obtain
Fj+1/2= Fj+1/2h + 12− W ¯Rj+1/2∆x sj− ˙ϕj, (3.11)
with Fj+1/2h defined as in (3.5b). Thus, the time derivative is included in the inhomogeneous flux. The resulting semidiscretisation reads
bW,jϕ˙j−1+ bC,jϕ˙j − aW,jϕj−1+ aC,jϕj− aE,jϕj+1= bW,jsj−1+ bC,jsj, (3.12)
with coefficients aW,jetc. defined in (3.7b). The extension of (3.8) is therefore an implicit ODE system
of the form
B ˙ϕ + Aϕ = Bs + b. (3.13)
To conclude, we have to apply a suitable time integration method to (3.13), for which we employ the ϑ-method [9].
4
Stability analysis
In this section we investigate stability of the ϑ-method applied to the implicit ODE system (3.13). For that purpose, we assume that K(x) = Const and ignore the source term and boundary conditions. An analysis of (3.13) for a constant coefficient advection-diffusion equation in terms of dissipation and dispersion is presented in [16].
For s = b = 0, the ϑ-method applied to (3.13) reads 1
∆tB ϕ
n+1− ϕn + (1 − ϑ)Aϕn+ ϑAϕn+1= 0, (4.1)
where 0 ≤ ϑ ≤ 1. Stability of (4.1) is estabished in the following theorem. Theorem.
The ϑ-method (4.1) is stable if the following conditions hold (
0 ≤ ϑ < 12 and d ≤ 2(1−2ϑ)1 C(R),
1
4 STABILITY ANALYSIS 8
with d = ∆t/∆x2, R = w∆x/K and where C is defined by C(z) = 1 − 1 − 2W (z)
1 − zW (z) 1 + z 12− W (z)2 , see Figure 5.
Proof. Written componentwise, the scheme (4.1) reads ∆x 12− W (R) 1 ∆t ϕ n+1 j−1 − ϕnj−1 + ∆x 12 + W (R) 1 ∆t ϕ n+1 j − ϕnj+ (1 − ϑ)K ∆x − B −ϕn j−1+ B−+ B+ϕnj − B+ϕnj+1+ ϑK ∆x − B −ϕn+1 j−1 + B −+ B+ϕn+1 j − B+ϕn+1j+1 = 0, (4.2)
with B± := B(±R). Substituting the discrete Fourier mode ϕnj = gneiκxj,
we obtain the following expressions for the amplification factor g, i.e., g = a1− (1 − ϑ)db1+ i(a2− (1 − ϑ)db2)
a1+ ϑdb1+ i(a2+ ϑdb2)
(4.3a) a1 = 12 − W (R) cos ψ +12 + W (R), a2 = − 12 − W (R) sin ψ (4.3b)
b1= 2 B++ B− sin2 12ψ, b2= R sin ψ, ψ = κ∆x. (4.3c)
From the stability requirement |g|2 ≤ 1 we can deduce a1b1+ a2b2+ (ϑ −12)d b21+ b22 ≥ 0,
or equivalently, in terms of the vectors a = (a1a2)Tand b = (b1b2)T,
a·b + ϑ −12d|b|2 ≥ 0. (∗)
From (∗) we conclude that if the scheme is stable for ϑ = 12, then it is certainly stable for12 < ϑ ≤ 1. So, consider first θ = 12, for which (∗) reduces to a·b ≥ 0. Expressing a and b in terms of 12ψ, substituting the resulting formulae and using the relation 12 B(z) + B(−z) = 1 + z 12 − W (z) we obtain the inequality
1 − sin2 12ψ 1 − 2W (R)(1 − RW (R) ≥ 0,
which should hold for 0 ≤ ψ < π. Clearly, a sufficient condition reads 1 − 2W (R) 1 − RW (R) ≤ 1,
which is identically satisfied, implying stability for 12 ≤ ϑ ≤ 1. Next, consider the case 0 ≤ ϑ < 1
2. In
this case, a sufficient condition for (∗) to hold is 1 − 1 − 2W (R)
1 − RW (R) − (2 − 4ϑ)d 1 + R 1
2 − W (R)
2 ≥ 0,
5 BOUNDARY TREATMENT 9 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 z C(z)
Figure 5: The function C.
5
Boundary treatment
In this section we outline the boundary conditions for the semidiscretisation (3.12).
At the surface, plankton nor nutrients can leave or enter the domain, therefore the corresponding fluxes must be 0, i.e.,
wP (0, t) − K(0, t)∂P ∂x(0, t) = 0, ∂Z ∂x(0, t) = ∂N ∂x(0, t) = 0, t > 0. (5.1) The bottom of the mixed layer is an artificial boundary where conditions are less obvious. So we need numerical boundary conditions. First, consider the semidiscretisation (3.12) for ϕ = P at the boundary point xM = 100, i.e.,
bW,MP˙M −1+ bC,MP˙M− aW,MPM −1+ aC,MPM− aE,MPM +1= bW,MsM −1+ bC,MsM, (5.2)
where the coefficients aC,M, bC,M etc. are defined in (3.7b). Note that this equation contains the
un-known PM +1in the virtual grid point xM +1. In order to eliminate PM +1we assume that there is no
dif-fusion in the virtual gridpoint, i.e., we set KM +1= 0. Taking into account the definitions of (weighted)
average in (3.6) we conclude that aE,M = 0, and consequently, PM +1disappears naturally from (5.2).
On the other hand, substituting KM +1= 0 in (3.11) we obtain the equivalent condition
FP,M +1/2 = wPM+ 12∆x sM − ˙ϕM, (5.3)
for the numerical flux at the virtual interface xM +1/2, where the subscript P is added to denote the flux
corresponding to the plankton concentration, thus the numerical flux is determined by the balance of advection and net production at x = xM +1/2. Next, for zooplankton, we derive a similar condition.
In the absence of advection the semidiscretisation (3.12) for ϕ = Z at x = xM reduces to the central
difference discretisation ∆x ˙ZM − 1 ∆x ¯K M −1/2ZM −1− ( ¯KM −1/2+ ¯KM +1/2)ZM + ¯KM +1/2ZM +1 = ∆xsM. (5.4)
In this case the condition KM +1= 0 is not sufficient to eliminate ZM +1and therefore we put ¯KM +1/2 =
0. This is equivalent to the zero-flux condition FZ,M +1/2 = 0. Finally, the concentration of nutrients
below the mixed layer is assumed to be constant, which means that the following boundary condition holds
6 NUMERICAL EXAMPLES 10
6
Numerical examples
We have computed numerical solutions in the mixed layer for a period of 4 years using 100 grid points and 6000 time steps. We used ϑ = 12 for the time integration of (3.13), which proved to be uncondi-tionally stable. The parameter values used in the simulations can be found in Table 1 and are based on [2, 4, 12, 14]. The behavior after some time, say 1 or 2 years, is independent of the initial conditions. We have chosen P (x, 0) = Z(x, 0) = 1 and N (x, 0) = 0.01x. This means that N varies from 0 at the surface to 1 at the bottom of the mixed layer.
First, we investigate the case without seasonal variability: µK = µp = 0. The solution is shown in
Figure 6(a). We see that it converges to the stationary solution, which is shown in Figure 7. The profiles are similar to the ones in [7]. The concentration of phytoplankton is maximal at a depth of around 20 meters and is also quite large at the surface and at 50 meters depth. The concentration of nutrients is small near the surface and starts to increase significantly from 40 meters downward. The main difference is a smaller concentration of zooplankton near the surface in [7], which is due to a sinking term for zooplankton.
Making the diffusion seasonally variable, µK = 0.5, we obtain the results in Figure 6(b). In this case
we see that P, Z and N slowly vary with the diffusion. This is expected, because higher diffusion means that nutrients diffuse towards the surface more rapidly. More nutrients results in more phytoplankton and more phytoplankton results in more zooplankton.
When the phytoplankton production rate is seasonally variable, µp = 0.4, we obtain the results in
Figure 6(c) and 6(d). Here we observe an important effect: a very steep phytoplankton bloom occurs in March. Levels are back to normal approximately one month after the bloom starts, so it lasts only a
Parameters Value H 100 m tmax 4 years Kmax 1.2773 × 10−3m2s−1 w 10−4m s−1 pmax 1.8 × 10−5s−1 ν 0.2 m−1 Lin 6 × 10−4mol photons m−2s−1 LH 2 × 10−5mol photons m−2s−1 `P 10−6s−1 `Z 5 × 10−7s−1 G 3 × 10−6s−1 NH 10−2mmol N2m−3 kiv 16.67 mmol N2 −1 m3 Pmin 1.2 × 10−2mmol N2m−3 γ 0.5 α 0.0175 r 0.3 µK 0, 0.5 µp 0, 0.4
6 NUMERICAL EXAMPLES 11
(a) µK= 0, µp= 0
(b) µK= 0.5, µp= 0
(c) µK= 0, µp= 0.4
(d) µK= 0.5, µp= 0.4
6 NUMERICAL EXAMPLES 12 0 20 40 60 80 100 0 0.05 0.1 0.15 x [m] C [mmol N2 m −3 ] P Z N
Figure 7: Stationary solution in case of no seasonal variability.
(a) (b)
Figure 8: Concentration of phytoplankton versus nutrients, and zooplankton versus phytoplankton, at a depth of 20 meters, for the case µK = 0.5 and µp = 0.4. In the NP figure, a darker shade of red
represents a higher concentration of zooplankton. In the PZ figure, darker shade of blue represents a higher concentration of nutrients.
very short period. During the phytoplankton bloom, also a zooplankton bloom starts. The zooplankton concentration does not increase as rapidly and as much as the phytoplankton concentration, but the zooplankton bloom does last longer. These phenomena correspond to observations in [17], [14] and [4]. The ways in which phytoplankton, zooplankton and nutrients interact, in case µK = 0.5 and µp =
counter-7 CONCLUDING REMARKS 13
clockwise direction. This corresponds to our expectations, because in both cases the species on the vertical axis (predators) feeds on the species on the horizontal axis (prey). This means that we get Lotka-Volterra type behaviour. Generally, when there is a lot of food, i.e., when the system is relatively far to the right in the figures, the concentration of predators will increase. This means that when the system is at the right, it will tend to move upward. Meanwhile, because the concentration of predators increases, the concentration of prey will decrease, which causes the system to move to the left. When the system is at the top left there are a lot of predators but not enough prey, and therefore the number of predators will decrease and the systems moves downward again. After enough predators have died, the concentration of prey can increase again, which means that the system moves to the right and the cycle is closed.
The phytoplankton bloom can easily be recognized in Figure 8. In Figure 8(a) it is the vertical peak and in Figure 8(b) it is the horizontal extreme. In both cases the bloom makes up a very large part of the curve. However, the trajectory moves much more rapidly in this part of the curve than in the rest of it. This must be the case, because the phytoplankton bloom lasts only about one month whereas the complete orbit corresponds to one year.
We see that two conditions are satisfied at the moment the phytoplankton bloom starts. First, the concentration of zooplankton is low and second, the concentration of nutrients is high. It is clear that the concentration of phytoplankton should increase in this case. However, the question remains why it increases so suddenly and why the bloom always occurs in March, regardless of the time of the year at which the simulation is started. The main reason appears to be that in March the water temperature increases and therefore the growth rate of the phytoplankton increases. The consequence of this increased growth rate is that the zooplankton, which is at that time only present in relatively small numbers, can not eat the phytoplankton quickly enough to keep the growth under control. This uncontrolled growth is approximately exponential, typical for an exploding population. The growth lasts until the concentration of nutrients becomes too small to support the growth of the phytoplankton and the concentration of zooplankton has become large enough.
7
Concluding remarks
The model that we propose describes the basic features of the ecosystem. One thing that makes this model different from most others, is the way we deal with the boundary at the bottom of the domain: we assume that the diffusion of phytoplankton and zooplankton is 0 at the lower boundary. This provides a good approximation of the real behavior. Another important aspect of the model is that we take into account seasonal variability. The seasonal effects of a higher phytoplankton production rate when the water is warmer and a higher turbulent diffusivity when the water is colder are included into the model.
The numerical scheme which we use to solve the system is the finite volume-complete flux scheme. The complete flux scheme includes the source term and time derivative already in the flux approximation and has proven to be a very accurate scheme.
We have applied the model to a (hypothetical) representative real life problem. Our results correspond to those in some of the literature. The results in other publications can be quite different, which is not surprising because (very) different parameter values are used. In order to validate the model properly, a geographical location should be chosen for which the parameter values can be obtained. The real behaviour at that location should be measured and the results of our model should be compared to that.
A few extensions of our model are suggested. It would be useful to have a separate model to calculate the hydrodynamical behaviour. This way it is possible to get more accurate profiles for, e.g., the tem-perature and the diffusion coefficient. Another extension would be to include transport in the horizontal
REFERENCES 14
direction. Doing this, advection will play a bigger role, but this should not be a problem considering the numerical scheme that is used. Another possibility is to include more species, such as more kinds of plankton.
References
[1] O. Arino, K. Boushaba and A. Boussouar (2000). A mathematical model of the dynamics of the phytoplankton-nutrient system. Nonlinear analysis: real world applications 1, 69-87.
[2] A. Dube and G. Jayaraman (2008). Mathematical modelling of the seasonal variability of plankton in a shallow lagoon. Nonlinear Analysis 69, 850-865.
[3] H.I. Freedman (1980). Deterministic Mathematical Models in Population Ecology, Marcel Dekker, New York.
[4] J.A. Freund, S. Mieruch, B. Scholze, K. Wiltshire and U. Feudel (2006). Bloom dynamics in a seasonally forced phytoplankton-zooplankton model: trigger mechanisms and timing effects. Ecological Complexity3, 129-139.
[5] L. Grieco, L.B. Tremblay and E. Zambianchi (2005). A hybrid approach to transport processes in the Gulf of Naples: an application to phytoplankton and zooplankton population dynamics. Continental Shelf Research25, 711-728.
[6] A. Larsson (2004). The role of mixing in an Antarctic ocean ecosystem: observations and model computations of vertical distributions of related parameters. Deep Sea Research Part II: Topical Studies in Oceanography51, 2807-2825.
[7] D.M. Lewis (2005). A simple model of plankton population dynamics coupled with a LES of the surface mixed layer. Journal of Theoretical Biology 234, 565-591.
[8] J.F. Lopes, A.C. Cardoso, M.T. Moita, A.C. Rocha and J.A. Ferreira (2009). Modelling the tem-perature and the phytoplankton distributions at the Aveiro near coastal zone, Portugal. Ecological Modelling220, 940-961.
[9] R.M.M. Mattheij, S.W. Rienstra and J.H.M. ten Thije Boonkkamp (2005). Partial Differential Equations, Modeling, Analysis, Computation. SIAM, Philadelphia.
[10] http://earthobservatory.nasa.gov/Features/Phytoplankton.
[11] S.V. Patankar (1980). Numerical Heat Transfer and Fluid Flow. Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere Publishing Corporation, New York.
[12] N.N. Pham Thi (2005). On positive solutions in a phytoplankton-nutrient model. Journal of Com-putational and Applied Mathematics117, 467-473.
[13] N.N. Pham Thi, J. Huisman and B.P. Sommeijer (2005). Simulation of three-dimensional phy-toplankton dynamics: competition in light-limited environments. Journal of Computational and Applied Mathematics174, 57-77.
REFERENCES 15
[14] N. Skliris, K. Elkalay, A. Goffart, C. Frangoulis and J.H. Hecq (2001). One-dimensional modelling of the plankton ecosystem of the north-western Corsican coastal area in relation to meteorological constraints. Journal of Marine Systems 27, 337-362.
[15] J.H.M. ten Thije Boonkkamp and M.J.H. Anthonissen (2011). The finite volume-complete flux scheme for advection-difussion-reaction equations. J Sci Comput 46, 47-70.
[16] J.H.M. ten Thije Boonkkamp and M.J.H. Anthonissen (2010). Extension of the complete flux scheme to time-dependent conservation laws. In: G. Kreiss et al. (eds.) Numerical Mathematics and Advanced Applications 2009, Proceedings ENUMATH 2009, Springer-Verlag, Berlin, 865-873.
[17] C. Troupin (2010). Seasonal variability of the oceanic upper layer and its modulation of biological cycles in the Canary Island region. Journal of Marine Systems 80, 172-183.
[18] A. Verma, S. Goorden, E. Alemayehu and P.A. Marin Zapata (2010). Modeling of the residence time of phytoplankton in the surface mixed layer of the ocean. private communication.
[19] J. Woods (2005). Stability and predictability of a virtual plankton ecosystem created with an individual-based model. Progress in Oceanography 67, 43-83.
