CWI  
  =-

Three-dimensional Simulation of Phytoplankton Dynamics

Exacte Wetenschappen
 
Programma Computational Science

1. General Project Information

1. 1. Project Title

The project is entitled: Three-dimensional Simulation of Phytoplankton Dynamics

1. 2. Project Acronym

The title is sufficiently short to avoid the use of an acronym.

1. 3. Principal Investigator

Dr. B.P. Sommeijer
Centre for Mathematics and Computer Science (CWI)
Department Modelling, Analysis and Simulation
Kruislaan 413, 1098 SJ Amsterdam
P.O.Box 94079, 1090 GB Amsterdam
tel. 020 - 592 4192, fax 020 - 592 4199
e-mail: bsom@cwi.nl

In collaboration with Dr. J. Huisman, Aquatic Microbiology, UvA

Dr. B.P. Sommeijer will act as contact person.

2. Summary

Phytoplankton takes up CO2 from the atmosphere and provides the basis of the food chains in the oceans. Modelling of phytoplankton is vital to understanding the dynamics of phytoplankton populations and the feedback of marine phytoplankton on the greenhouse effect.

For many years, phytoplankton models have been studied in a well-mixed environment, leading to the common paradigm that phytoplankton blooms cannot develop if the upper water layer of a stratified water column is deeper than some critical depth [25]. Huisman et al. [16] discovered that also the vertical mixing rate is a critical parameter for the development of blooms.
Using a diffusion-reaction model, they found that if the vertical diffusion is less than a critical value, a bloom can develop, irrespective of the depth of the water column. Recently, these one-dimensional models were extended by including vertical transport processes [5, 14]. This is motivated by the fact that phytoplankton species may have a specific weight larger/smaller than water, causing sinking/buoyancy of these species. These studies revealed many interesting and surprising results.

In all the above models, only variation in the vertical direction was considered. An important aim of this project is to incorporate horizontal flow of the water. The resulting three-dimensional PDEs will be used to study the so-called `biological pump'. By this we mean the sequestration of CO2 from the atmosphere which is then transported to the bottom of the ocean by sinking phytoplankton species. This mechanism significantly contributes to the reduction of the greenhouse effect. First results obtained by the 1D model provide new insight in the relevant parameters [18]. Furthermore, we will include a module for the nutrients (which, up to now, were assumed to be amply available), and we will take into account the influence of water column stratification on phytoplankton dynamics.
The combined model is a system of 3D integro-partial differential equations, requiring novel and tailor-made numerical algorithms.

3. Classification

The research topic can be classified in the following disciplines:
  1. Numerical mathematics;
  2. Population dynamics.

4. Composition of the Research Team

The research team will be composed as follows:
ResearcherSpecializationFinanced by Fte/year
dr. B.P. Sommeijernumerical analysisCWI0.3
dr. J. Huismanbiology/modelUvA0.1
prof. dr. J.G. Verwernumerical analysisCWIpm.
drs. J. Koksoftware/HPC/visualizationCWI0.1
drs. J. Passargebiology/modelUvA0.1
OIO, NNalgorithms/softwareEW, Comp. Sc.1.0

We request funding for a Ph.D. student (OIO). He/she should have a strong background in numerical methods for partial differential equations and ordinary differential equations and, preferably, in parallel computing. Also some expertise in biological modelling is very welcome.

The project should result in a Ph.D. thesis. As promotors will act prof. dr. J.G. Verwer (for the numerical aspects) and dr. J. Huisman (for the biological/modelling aspects). It is to be expected that dr. Huisman will be shortly appointed as professor of Aquatic Microbiology at the University of Amsterdam (UvA). Prof. Verwer is head of the department Modelling, Analysis and Simulation at CWI and professor of numerical analysis at the University of Amsterdam. Dr. Huisman is the new leader of the department of Aquatic Microbiology at the Institute for Biodiversity and Ecosystem Dynamics of the UvA and expert in the field of phytoplankton dynamics.
The integration of the specific expertise of numerical mathematicians (from the CWI-MAS group) and biologists (from the UvA-IBED group) can lead to a major step forward in a better understanding of the dynamics of aquatic ecosystems. Therefore, a close co-operation is of paramount relevance for the success of the project.

5. Research School

The members of the research team actively participate in two research schools: Thomas Stieltjes School for Mathematics, for MAS-CWI, and Research School SENSE (Socio-Economic and Natural Sciences of the Environment), for Aquatic Microbiology, UvA.

6. Description of the Proposed Research

6. 1. Scientific relevance of phytoplankton models

Modelling the dynamics of phytoplankton is of great importance to many aspects of human interest. The most important aspect is the role that sinking phytoplankton species play in the downward export of biological carbon to the bottom of the ocean [24, 23], the so-called `biological pump'. By this mechanism, several gigatons per year of carbon dioxide are removed from the atmosphere [7], thus making a significant contribution to the reduction of the greenhouse problem on earth. In fact, it has been suggested to increase the phytoplankton concentration in the upper water layer by nutrient enrichment of the ocean. In this way it is hoped to obtain an enhanced sequestration of carbon from the atmosphere. Some first experiments in this direction, by means of iron-fertilization of the Southern Ocean, have taken place. These experiments indicate, however, that the measured amount of downward carbon export hardly increases, in spite of the enhanced phytoplankton concentration in the upper layer [3]. Based on predictions resulting from our model analysis, we describe in [18] a possible explanation for this failure. It is made plausible that the increased turbidity of the water plays a crucial role. However, the models that we are using are still quite limited and require significant refinement in order to produce predictions that are sufficiently reliable for policy makers to base decisions on.

6. 2. Need for supercomputing and the role of numerical mathematics

Along with the two classical scientific lines of experimentation and theoretical analysis, the value of computer simulation of phytoplankton dynamics has now been established. The ever-increasing capacity of computers has given an enormous impulse to the development of codes for performing real-life simulations. The complexity of the full system, however, is so large, that a number of simplifications still has to be incorporated in the current models, simply to keep the amount of storage and CPU at a realistic level. Nevertheless, a possibility to bridge part of the gap between present-day practice and the ultimate goal is to use efficient, tailor-made numerical algorithms and innovative computer science techniques. In this way, satisfactory results are feasible to make a significant step towards understanding the full complexity of the system.

Initially, the development and testing of the algorithms will be performed on workstations (which are available at CWI). At a later stage the simulation of real-life situations will be carried out on a parallel machine, like the new SGI Origin 3000 (`Teras') at SARA.

6. 3. Problem Formulation and Anticipated Results

6. 3. 1. Prior research

Critical parameters. In recent years, research on phytoplankton models has received much attention. Initially, only models in a well-mixed environment were studied. This led to the conclusion that phytoplankton blooms can develop only in waters where the depth of water-column stratification is shallower than some critical depth (see e.g. [25]). This can be understood by noticing that only in shallow water the light conditions are sufficient to let (depth-integrated) photosynthesis exceed the (depth-integrated) natural losses of phytoplankton. However, field observations did not always confirm this critical depth concept [27, 6], but a firm theory to explain these observations was lacking. A crucial step forward was to release the assumption of vertical homogeneity and to study the phytoplankton distribution over depth [16]. As a result of this approach, the model changes to an integro-partial differential equation, with variation in time and in space. A surprising result of this study is that there exist two different mechanisms that may cause a phytoplankton bloom. Apart from the critical depth, it was found that the diffusive mixing rate plays a crucial role: below a certain critical mixing rate (i.e., in quiet waters) blooms can always develop, independent of the water depth.
Transport phenomena. A next step is to include vertical transport processes (other than by diffusion) into the models. This is motivated by the fact that phytoplankton species often have a specific weight different from that of water. Species with a specific weight larger than water will sink downwards. On the other hand, species with a smaller weight give rise to buoyancy. Since phytoplankton requires light for its photosynthesis, buoyant species have an advantage over sinking species, because buoyant species can easily reside in the well-lit zone (i.e., close to the water surface). Therefore, including the vertical velocity of specific species in the model has a significant impact. In mathematical terms we have now changed to a advection-diffusion-reaction model. First results have been reported in [5, 14, 19].
Competition models. In each of the aforementioned studies only one particular phytoplankton species has been considered. A natural step is to study systems in which various species are considered in a competition context. An important factor that will play a role is that species that stay in the upper layer (and hence, capture a lot of light) have the additional advantage to shade other species at deeper levels. Nevertheless, as it turns out, the balance between survival and depth for the various species is very delicate and depends on values of other parameters as well, such as background turbidity, incident light intensity, and specific light attenuation coefficients of the different species [17].

Competition models have been studied before by many authors, but always using simpler models (see [26, 9, 15]). Recent work on horizontal turbulence in a competition context seems very promising [4, 20].

6. 3. 2. Proposed research

All topics described in the preceding subsection are very promising, in the sense that new, and sometimes surprising, results were obtained with the model. However, a common restriction in all applications so far is that either only variation in the vertical direction has been considered or only variation in the horizontal direction. The main aim of the current project is to combine the effects due to the various spatial directions. Hence, we will study the full three-dimensional growth and transport of phytoplankton, for example by imposing a horizontal wind field or taking into account differences in barometric pressure.

A priori, the influence of such effects is hard to predict, but it will certainly play a significant role. For example, it is already well-known that horizontal turbulence may disrupt bloom formation if horizontal dispersal rates exceed local growth rates (see e.g., [21]). Little is known, however, of the interactive effects of horizontal and vertical turbulent diffusion on the dynamics of sinking or buoyant phytoplankton, i.e. the three-dimensional structure of phytoplankton blooms. Generally speaking, simple flows in three dimensions can potentially generate very complicated patterns [8].

Another important feature for the dynamics is the influence of nutrients. In all above studies it was tacitly assumed that an ample supply of nutrients was available. It is of course more realistic to release this ideal nutrient condition and to study combined nutrient-phytoplankton models. Such models have been studied before (see e.g., [26, 9]), but usually in a well-mixed, homogeneous environment. Hence, the spatial variation that is incorporated in our models gives (literally) `an extra dimension'. This topic nicely fits in the 3-dimensional model extension proposed above since nutrient transport will be largely determined by the horizontal flow as well as by upwelling from below.

A third relevant issue, as already pointed out in Section 6.1, is the role of the turbidity of the water column. In [18] it was shown that the maximal sinking velocity that can be sustained by a sinking phytoplankton population is inversely proportional to the background turbidity. Consequently, an increased turbidity (for example by iron-fertilization) will lead to a decreased maximal vertical velocity, which may have a negative influence on the export of carbon to the deep ocean. Since this export is directly related to the greenhouse problem, it is of high relevance for society.

Finally, we want to develop models to study the competition between many phytoplankton species. In [17], some results have been published based on a model that did not take into account the effects of sinking and buoyancy. Moreover, only one spatial dimension (i.e., the vertical) was considered and only a few resources. This model predicts a low species diversity: blooms are dominated by one or a few species only. It is highly interesting to study whether this prediction remains valid in the more sophisticated models where three spatial dimensions and vertical transport have been incorporated, and where both the number of species and the number of resources will be increased.

To summarize: the new modelling issues of this proposal are:

  1. the combination of horizontal and vertical variability;
  2. a combined three-dimensional nutrient-phytoplankton model;
  3. the role of turbidity on carbon export by sinking phytoplankton;
  4. three-dimensional competition models.
Numerical algorithms. As said before, the mathematical model is of advection-diffusion-reaction type with a highly nonlinear `reaction' term originating from the space-dependent light gradient. This reaction term, which is in fact the production term generated by phytoplankton photosynthesis, appears in the model as an integral (over depth) of the unknown concentrations. Hence, the model is of so-called `integro-partial differential equation of Fredholm type' in 3D and for its solution in the space-time domain we need a numerical approach. We will follow the Method of Lines approach, that is, first the spatial differential operators, as well as the integral term, will be replaced by discrete approximations and subsequently the huge resulting system of ordinary differential equations (ODEs) -- which is still continuous in time -- will be integrated numerically.

For the spatial discretization we will use a (conservative) finite volume technique based on so-called upwind approximations for the advection terms and central approximations for the diffusion terms (see e.g., [11] for details). The integral in the reaction term will be approximated by a suitable quadrature formula. The population density over space is often far from uniform. This means that efficient use can be made of a non-uniform spatial grid: in regions with a high spatial activity a fine mesh will be used and a coarse mesh elsewhere. In this way the total number of grid points, and hence the CPU time, can be drastically reduced.

Next, the resulting system of ODEs is integrated in time. Because of the wide range in temporal solution scales, this ODE system is stiff and, to avoid a severe timestep restriction for stability reasons, an implicit integration technique seems to be an appropriate choice. In the literature (see e.g. [10]) a wide variety of implicit integration techniques has been described. However, in using this type of methods we are faced with the task to solve, in each time step, a system of implicit relations to obtain the solution at the next point in time. Due to the three spatial dimensions in combination with the integral term, the resulting linear algebra problem is extremely expensive. Therefore, standard implicit techniques such as Crank-Nicolson or BDF-type methods are not feasible and we have to change to a time integration process that is adapted to this particular application.

A possibility is to consider special time-stepping techniques, such as

A common feature of these methods is that they circumvent the coupled three-dimensional computations. In fact, this huge problem is replaced by a sequence of simpler ones. However, all these techniques are rather `special purpose' in nature and a successful use requires adaptation to the particular application. In our case of phytoplankton dynamics a thorough analysis should reveal whether the above techniques lead to a feasible solution process and, if so, which is the most efficient one. If none of these techniques is satisfactory, then new algorithms have to be developed.

Another feature which renders these models numerically interesting is the coupling of the dynamical model with the equations for the nutrients and for dead organic matter.
It may turn out that here an operator-splitting technique (of Strang type) will be most appropriate.

6. 4. Comparison with other groups

To the best of our knowledge, the approach to model the dynamics of phytoplankton by means of three-dimensional advection-diffusion-reaction equations, in the form of integro-PDEs is new and has not been done before. However, our own recent work on one-dimensional vertical plankton models [14, 16, 18] and recent work by others on two-dimensional horizontal plankton models [4, 20] suggest that a three-dimensional approach will be very promising, and should be given high priority if we wish to keep ahead in this exciting field of research.

6. 5. Co-operation with other groups

Inside the Netherlands Outside the Netherlands

6. 6. Embedding

In the department MAS at CWI, `Applications to the Life Sciences' is one of the major research topics. This new topic (started January 2000) focuses on mathematical problems from biology and medicine. Co-operation has been started with groups working in cell-, neuro- and microbiology. All models have been formulated in terms of PDEs, the numerical solution of which has been a major research topic at CWI for the last 25 years. Furthermore, during the last 10 years, expertise has been acquired in constructing parallel methods for PDEs.
Hence, the proposed research nicely fits in the field of interest of the department MAS.

For the department of Aquatic Microbiology, which is embedded in the Institute for Biodiversity and Ecosystem Dynamics (IBED) of the University of Amsterdam, the study of the dynamics of phytoplankton species is the main activity. This has been the case for many years, and certainly will be continued in the years to come since dr. Huisman recently received a `NWO-PIONIER' grant to further stimulate the research in this direction. The Aquatic Microbiology group combines mathematical modelling, laboratory experiments and field research on phytoplankton, spanning the entire range from molecular biology and physiology to the population dynamics and ecosystem impacts of phytoplankton. Hence, it is obvious that the proposed project is a natural part of the current activities of this group.

6. 7. Utilization aspects

7. Work Programme

The 4-year research programme comprises six activities:
  1. The first 6 months of the project are reserved for the OIO to get acquainted with the relevant literature on this subject, including the recent papers by Huisman et al.
  2. The second step will be to extend the one-dimensional model to a three-dimensional model. This step requires the construction and analysis of a suitable numerical solver, as well as its implementation. Since this extension is an important task in the project, it is expected to take a full year.
  3. Then we will investigate the influence of nutrients, when incorporated in this three-dimensional model. In this way, our usual assumption of an ample supply of nutrients is made more realistic. This activity is planned to require 6 months.
  4. Next, in a stratified version of the model, carbon export will be studied. Special attention will be paid to the role of the turbidity since it turned out to be a critical parameter in the one-dimensional model. This part is estimated to be completed within 9 months.
  5. The last topic concerns competition models. Here the models will be modified as to allow the study of multi-species dynamics. It is to be expected that the computing time and memory requirements will drastically increase in these kind of simulations. Therefore, the use of a parallel computer like the `Teras' is necessary in this part of the project. Also the numerical algorithms will probably need some modifications. This topic will take another year.
  6. The remaining time will be reserved for documentation and preparation of the thesis.

8. Expected Use of Instrumentation

  1. This project does not require special equipment.
  2. In the second part of the project the use of a supercomputer will be necessary (see also point 5 of Section 7). A total of 250 hours CPU time is a good estimate.

9. Requested Budget

The costs for the full 4-year period are specified in the following table.

PersonnelFte/yearEW, Comp. Sc. costs
 fl
dr. B.P. Sommeijer0.3 0
dr. J. Huisman0.1 0
prof. dr. J. G. Verwerpm. 0
drs. J. Kok0.1 0
drs. J. Passarge0.1 0
OIO, NN1.0 223,000
Travel budget  7,350
Consumables (software, etc.)  10,000
 
Totals  240,350

10. Literature

Five most relevant, recent publications of the research team
  1. P.J. van der Houwen and B.P. Sommeijer, Approximate factorization for time-dependent partial differential equations, J. Comput. Appl. Math. 128, 447 - 466 (2001).
  2. J.G. Verwer and B.P. Sommeijer, A numerical study of mixed parabolic-gradient systems, J. Comput. Appl. Math. 132, 191 - 210 (2001).
  3. J. Huisman and B.P. Sommeijer, Population dynamics and export production of sinking phytoplankton in stratified waters, final version in preparation, CWI Report, MAS-R01xx, submitted to Limnol. Oceanogr. (2001).
  4. J. Huisman, P. van Oostveen and F.J. Weissing, Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms, Limnology and Oceanography 44, 1781-1788 (1999).
  5. J. Huisman and F.J. Weissing, Biodiversity of plankton by species oscillations and chaos, Nature 402, 407-410 (1999).

References

  1. U.M. Ascher, S.J. Ruuth and B. Wetton, Implicit-explicit methods for time-dependent PDEs, SIAM J. Numer. Anal. 32, 797-823 (1997).
  2. H.J.W. de Baar, et al., Importance of iron for plankton blooms and carbon-dioxide drawdown in the Southern Ocean, Nature 373, 412-415 (1995).
  3. P.W. Boyd, et al., A mesoscale phytoplankton bloom in the polar Southern Ocean stimulated by iron fertilization, Nature 407, 695-702 (2000).
  4. A. Bracco, A. Provenzale and I. Scheuring, Mesoscale vortices and the paradox of the plankton, Proc. R. Soc. Lond. B 267, 1795-1800 (2000).
  5. U. Ebert, M. Arrayás, N.M. Temme, B.P. Sommeijer and J. Huisman, Critical conditions for phytoplankton bloom, final version in preparation, CWI Report, MAS-R01xx, submitted to J. of Math. Biology (2001).
  6. H.C. Eilertsen, Spring blooms and stratification, Nature 363, 24 (1993).
  7. P.G. Falkowski, R.T. Barber and V. Smetacek, Biogeochemical controls and feedbacks on ocean primary production, Science 281, 200-206 (1998).
  8. G.O. Fountain, D.V. Khakhar, I. Mezic and J.M. Ottino, Chaotic mixing in a bounded three-dimensional flow, J. Fluid. Mech. 417, 265-301 (2000).
  9. J.P. Grover, Resource Competition, Chapman & Hall, London (1997).
  10. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II; Stiff and Differential-algebraic Problems, Springer Series in Computational Mathematics 14, Springer-Verlag, Berlin, 2nd edition (1996).
  11. C. Hirsch, Numerical Computation of Internal and External Flows, I: Fundamentals of Numerical Discretization, Wiley and Sons, Chichester (1988).
  12. P.J. van der Houwen and B.P. Sommeijer, Splitting methods for three-dimensional transport models with interaction terms, J. Scient. Comput. 12, 215-231 (1997).
  13. P.J. van der Houwen and B.P. Sommeijer, Approximate factorization for time-dependent partial differential equations, J. Comput. Appl. Math. 128, 447 - 466 (2001).
  14. J. Huisman, M. Arrayás, U. Ebert and B.P. Sommeijer, How do sinking phytoplankton species manage to persist in unstratified waters?, CWI Report, MAS-R0105, submitted to Am. Nat. (2001).
  15. J. Huisman, R.R. Jonker, C. Zonneveld and F.J. Weissing, Competition for light between phytoplankton species: experimental tests of mechanistic theory, Ecology 80, 211-222 (1999).
  16. J. Huisman, P. van Oostveen and F.J. Weissing, Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms, Limnol. Oceanogr. 44, 1781-1788 (1999).
  17. J. Huisman, P. van Oostveen and F.J. Weissing, Species dynamics in phytoplankton blooms: incomplete mixing and competition for light, Am. Nat. 154, 46-68 (1999).
  18. J. Huisman and B.P. Sommeijer, Carbon export by marine phytoplankton requires clear water, Internal Report.
  19. J. Huisman and B.P. Sommeijer, Population dynamics and export production of sinking phytoplankton in stratified waters, final version in preparation, CWI Report, MAS-R01xx, submitted to Limnol. Oceanogr. (2001).
  20. G. Karolyi, A. Pentek, I. Scheuring, T. Tel and Z. Toroczkai, Chaotic flow: the physics of species coexistence, Proc. Natl. Acad. Sci. 97, 13661-13665 (2000).
  21. H. Kierstead and L. Slobodkin, The size of water masses containing plankton blooms, J. Mar. Res. 12, 141-147 (1953).
  22. G.I. Marchuk, Splitting and alternating direction methods, In: Handbook of Numerical Analysis I (Eds: P.G. Ciarlet and J.L. Lions), North-Holland, Amsterdam (1990).
  23. J.L. Sarmiento, T.M.C. Hughes, R.J. Stouffer and S. Manabe, Simulated response of the ocean carbon cycle to anthropogenic climate warming, Nature 393, 245-249 (1998).
  24. J.L. Sarmiento and C. le Quéré, Ocean carbon dioxide uptake in a model of century-scale global warming, Science 274, 1346-1350 (1996).
  25. H.U. Sverdrup, On conditions for the vernal blooming of phytoplankton, J. Cons. Perm. Int. Explor. Mer 18, 287-295 (1953).
  26. D. Tilman, Resource Competition and Community Structure, Princeton Univ. Press, Princeton (1982).
  27. D.W. Townsend, M.D. Keller, M.E. Sieracki and S.G. Ackleson, Spring phytoplankton blooms in the absence of vertical water column stratification, Nature 360, 59-62 (1992).
  28. J.G. Verwer, W.H. Hundsdorfer and J.G. Blom, Numerical time integration for air pollution models, to appear in Surveys on Math. for Industry (2001).

=-

Document http://www.cwi.nl , last modified: Jul 10 11:58