Bfm V5.1.0 Manual R1.1 201508

User Manual:

Open the PDF directly: View PDF .
Page Count: 104

Download
Open PDF In Browser	View PDF

The Biogeochemical Flux Model (BFM)
Equation Description and User Manual
BFM version 5.1

The BFM System Team

Release 1.1, August 2015
—– BFM Report series N. 1 —–
http://bfm-community.eu
info@bfm-community.eu

This document should be cited as:
Vichi M., Lovato T., Lazzari P., Cossarini G., Gutierrez Mlot E., Mattia G., Masina S., McKiver W. J., Pinardi N., Solidoro C., Tedesco L., Zavatarelli M. (2015). The Biogeochemical
Flux Model (BFM): Equation Description and User Manual. BFM version 5.1. BFM Report
series N. 1, Release 1.1, August 2015, Bologna, Italy, http://bfm-community.eu, pp. 104

Copyright 2015, The BFM System Team. This work is licensed under the Creative Commons
Attribution-Noncommercial-No Derivative Works 2.5 License. To view a copy of this license, visit
http://creativecommons.org/licenses/by-nc-nd/2.5/ or send a letter to Creative Commons, 171 Second Street, Suite 300,
San Francisco, California, 94105, USA.

2

Contents
I.

The BFM Equations

13

1. The formalism of the BFM equations

15

2. The pelagic plankton model
2.1. The environmental parameters affecting biological rates . . . . . . .
2.2. Phytoplankton . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1. Photosynthesis and carbon dynamics . . . . . . . . . . . . .
2.2.2. Multiple nutrient limitation . . . . . . . . . . . . . . . . . .
2.2.3. Nutrient uptake . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3.1. Nitrogen . . . . . . . . . . . . . . . . . . . . . .
2.2.3.2. Phosphorus . . . . . . . . . . . . . . . . . . . .
2.2.3.3. Silicate . . . . . . . . . . . . . . . . . . . . . . .
2.2.4. Nutrient loss associated with lysis . . . . . . . . . . . . . .
2.2.5. Exudation of carbohydrates . . . . . . . . . . . . . . . . .
2.2.6. Chlorophyll synthesis and photoacclimation . . . . . . . . .
2.2.7. Light limitation and photoacclimation based on optimal light
2.2.8. Iron in phytoplankton . . . . . . . . . . . . . . . . . . . . .
2.2.9. Phytoplankton sinking . . . . . . . . . . . . . . . . . . . .
2.3. Bacterioplankton . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1. Regulating factors . . . . . . . . . . . . . . . . . . . . . .
2.3.2. Respiration . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3. Mortality . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.4. BACT1 parameterization . . . . . . . . . . . . . . . . . .
2.3.4.1. Substrate uptake . . . . . . . . . . . . . . . . . .
2.3.4.2. Nutrient release and uptake . . . . . . . . . . . .
2.3.4.3. Excretion . . . . . . . . . . . . . . . . . . . . . .
2.3.5. BACT2 parameterization . . . . . . . . . . . . . . . . . . .
2.3.5.1. Substrate uptake . . . . . . . . . . . . . . . . . .
2.3.5.2. Nutrient release and uptake . . . . . . . . . . . .
2.3.6. BACT3 parameterization . . . . . . . . . . . . . . . . . .
2.3.6.1. Substrate uptake . . . . . . . . . . . . . . . . . .
2.3.6.2. Nutrient release and uptake . . . . . . . . . . . .
2.3.6.3. Excretion . . . . . . . . . . . . . . . . . . . . . .
2.4. Zooplankton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1. Food availability . . . . . . . . . . . . . . . . . . . . . . .
2.4.2. Ingestion . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3. Excretion/egestion . . . . . . . . . . . . . . . . . . . . . .

19
19
20
22
24
25
25
25
26
26
26
27
28
29
29
29
31
31
31
33
33
33
34
34
34
34
34
34
35
35
35
35
37
37

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

3

Contents
2.4.4. Respiration . . . . . . . . . . . . . . .
2.4.5. Excretion/egestion of organic nutrients
2.4.6. Inorganic nutrients . . . . . . . . . . .
2.5. Non-living components . . . . . . . . . . . . .
2.5.1. Oxygen and anoxic processes . . . . .
2.5.2. Dissolved inorganic nutrients . . . . . .
2.5.3. Dissolved and particulate organic matter
2.5.4. The carbonate system . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

38
39
39
40
40
41
43
44

3. The sea ice biogeochemical model
3.1. A model for sea ice biogeochemistry . . . . . . . . . . . .
3.2. Model structure . . . . . . . . . . . . . . . . . . . . . . .
3.3. Sea ice Algae Dynamics . . . . . . . . . . . . . . . . . .
3.4. Nutrient Supply and Dynamics . . . . . . . . . . . . . . .
3.5. Gases and Detritus . . . . . . . . . . . . . . . . . . . . .
3.6. The coupling strategy . . . . . . . . . . . . . . . . . . . .
3.6.1. Boundary fluxes: the sea ice-ocean interface . . . .
3.6.2. Boundary fluxes: the sea ice-atmosphere interface

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

49
49
49
50
55
56
56
57
59

II. BFM code description
4. Installation, configuration and compilation
4.1. Installation . . . . . . . . . . . . . . . . .
4.1.1. System Requirements . . . . . . .
4.2. Configuration: BFM Presets . . . . . . . .
4.3. Compiling the BFM STANDALONE . . .

61

.
.
.
.

63
63
63
65
67

.
.
.
.
.
.
.
.
.
.
.
.
.

69
69
69
69
71
71
73
73
73
73
73
76
76
76

6. Structure of the code
6.1. Coding rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2. Memory layout and variable definition . . . . . . . . . . . . . . . . . . . . . . . . .

81
81
81

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

5. Running the STANDALONE model
5.1. Description . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2. State variables initial conditions (IC) . . . . . . . . . . . . . .
5.3. Numerical integration . . . . . . . . . . . . . . . . . . . . . .
5.4. Forcing functions . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1. Analytical forcing functions . . . . . . . . . . . . . .
5.4.2. Boundary forcing for atmospheric CO2 . . . . . . . .
5.4.3. Environmental data from file . . . . . . . . . . . . . .
5.4.4. External event data . . . . . . . . . . . . . . . . . . .
5.5. Test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5.1. STANDALONE_CO2TEST: carbonate system model
5.5.2. STANDALONE_PELAGIC: the full pelagic BFM . .
5.5.3. STANDALONE_SEAICE: sea ice biogeochemistry . .
5.6. Output visualization . . . . . . . . . . . . . . . . . . . . . . .

4

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

Contents
6.3. The BFM flow-chart . . . . . . . . . . . . .
6.3.1. Initialization procedures . . . . . . .
6.3.2. Time marching . . . . . . . . . . . .
6.3.3. Computation of pelagic reaction terms

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

83
83
84
87

7. Design model layout components
7.1. Layout syntax and modular structure . . . . . . . . . . . .
7.2. Example 1. Adding a subgroup . . . . . . . . . . . . . . .
7.3. Example 2. Adding a biochemical component and a group
7.4. Example 3. Removing components from zooplankton . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

89
89
92
94
95

8. Model outputs
8.1. Output . . . . . . . . . . . . . .
8.2. Restart . . . . . . . . . . . . . .
8.3. Aggregated diagnostics and rates
8.4. Mass conservation . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

97
97
97
97
99

Bibliography

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

101

5

List of Equation Boxes
2.1.
2.2.
2.3.
2.4.
2.5.
2.6.
2.7.

Phytoplankton . . . . . . . . . . . . . .
Chlorophyll synthesis . . . . . . . . . .
Bacteria . . . . . . . . . . . . . . . . .
Zooplankton . . . . . . . . . . . . . . .
Dissolved inorganic nutrient equations. .
Dissolved organic matter . . . . . . . .
Particulate organic detritus equations. .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

21
28
33
36
42
43
45

3.1. Sea ice algae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

7

List of Tables
1.1. BFM reference state variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. List of abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17
18

2.1.
2.2.
2.3.
2.4.
2.5.
2.6.
2.7.
2.8.
2.9.

.
.
.
.
.
.
.
.
.

21
23
30
32
36
37
39
42
47

3.1. List of the sea ice model state variables . . . . . . . . . . . . . . . . . . . . . . . .
3.2. Ecological and physiological parameters in BFM-SI. . . . . . . . . . . . . . . . . .

52
53

4.1. List of optional arguments for bfm_configure . . . . . . . . . . . . . . . . . . . . .

68

5.1.
5.2.
5.3.
5.4.
5.5.
5.6.

.
.
.
.
.
.

70
70
71
72
72
74

7.1. Example of model layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

Environmental variables . . . . . . . . . . .
Phytoplankton parameters . . . . . . . . . .
Iron cycle parameters . . . . . . . . . . . . .
Bacterioplankton . . . . . . . . . . . . . . .
Microzooplankton . . . . . . . . . . . . . . .
Mesozooplankton . . . . . . . . . . . . . . .
Switches for mesozooplankton food limitation
Chemical parameters . . . . . . . . . . . . .
Carbonate system . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

namelist bfm_nml in BFM_General.nml . . . . . . . .
namelist Param_parameters in BFM_General.nml . . .
namelist PAR_parameters in Pelagic_Environment.nml
namelist standalone_nml in Standalone.nml . . . . . .
namelist time_nml in Standalone.nml . . . . . . . . .
namelist forcings_nml in Standalone.nml . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

9

Introduction
Synopsis
This report is divided in two main parts. The first
part contains an extended description of the Biogeochemical Flux Model (BFM version 5) equations and as such it describes different parameterizations that can be used to simulate biogeochemical processes in the marine system. Part II
describes the model set-up, compilation and execution of the standard test cases. It also provides
the major technical implementation of the model
and information to understand the code structure.
This part gives examples of the possible diagnostics that can be activated and illustrates advanced
features that allow the user to modify the model
adding new state variables and components.
The BFM is a direct descendant of the European Regional Seas Ecosystem Model (ERSEM)
and shares most of its characteristics with
the original formulations (Baretta et al., 1995;
Baretta-Bekker et al., 1997). The major difference is that BFM focuses more on the biogeochemistry of lower trophic levels in marine
ecosystems than on trophic food webs, which
is where it takes its name. The philosophy
of the BFM and the mathematical formulation
applied throughout the book are described in
Vichi et al. (2007b). The users familiar with
ERSEM will find that the code notation described by Blackford and Radford (1995) has not
changed substantially and that the state variables
are indicated with the same naming convention.
Vichi et al. (2007b) have extended this notation
making it formally consistent and additional variables have been named accordingly.
The structure of the code was instead extensively modified with respect to the original
ERSEM code. The major reason for a full rewriting lied in the growing need to couple biogeochemical processes with hydrodynamical model
of various forms and to allow a modular expan-

sion of the number and type of state variables in
modern coding standards. The necessity to have
a flexible system which can be embedded in the
existing state-of-the-art ocean general circulation
models (OGCM) required a re-organization of the
structure, yet the fundamentals of the coupling
strategy are very similar to the first coupled implementations. From BFM version 5.1 the coupled configurations are publicly supported. The
BFM module can be easily coupled with different ocean models and an example of the coupled configuration with the NEMO ocean model
(http://nemo-ocean.eu) is given in the companying manual of this release Vichi et al. (2015).

The BFM module
This report documents the BFM equations for the
standalone pelagic configuration of the code that
solves the reaction part of the general equations
for biogeochemical tracers in the marine environment. The state variables presented in Part I
are historically derived from the pelagic component of ERSEM and represent the standard set
of transformations for the major chemical constituents. The BFM by construction solves the
cycles of C, O, N, P, Si in the lower trophic levels
of marine ecosystems and it also allows the inclusion of the Fe cycle and carbonate system chemistry by means of user-controlled flags. Different parameterizations are proposed for the various
functional groups as they have been demonstrated
valid in certain model applications. Given the inherent empirical nature of biogeochemical modelling, the proposed parameterizations are given “as
they are” and it is up to the user to decide which
one to use given the specific application. The
modular implementation of the BFM allows to
implement new parameterizations and new state
variables.

11

List of Tables
The equations of Part I consider no vertical or
horizontal processes and the system is solved as a
set of ordinary differential equations. At the heart
of the BFM module is the possibility to increase
the number of components, as for instance implementing any number of phytoplankton groups
that share the same dynamical terms but differ in
term of parameter and parameterization choices.
The user interested in this feature, in the modification of the model layout and in the definition of
specific diagnostics will find in Part II all the required information. The reference test cases will
be presented in the STANDALONE implementation, which is a time-dependent box model with a
given depth where pelagic processes are assumed
homogeneous.

distribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of
the License, or (at your option) any
later version. The BFM is distributed
in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE. See the
GNU General Public License for more
details.You should have received a copy
of the GNU General Public License
along with the BFM core. If not, see
.

Licensing

Acknowledgments

The BFM is free software and is protected by the GNU Public License
http://www.gnu.org/licenses/gpl.html.
All
files of the BFM contain the following statement:

The code support and development is currently
done by the BFM System Team whose members are listed on the BFM web site (http://bfmcommunity.eu). The initial core of the BFM system software was developed by Piet Ruardij and
Marcello Vichi. The BFM system team is grateful
to Job Baretta, Hanneke Baretta-Bekker, Piet Ruardij and Wolfgang Ebenhoeh for their contribution to the science of biogeochemical modeling.

Copyright 2013, The BFM System
Team (info@bfm-community.eu)

This file is part of the BFM.
The BFM is free software: you can re-

12

Part I.

The BFM Equations

13

1. The formalism of the BFM equations
The BFM equations are written following the
notation proposed by Vichi et al. (2007b), which
is briefly summarized in this section. The reader
is referred to Vichi et al. (2007b) for a more theoretical viewpoint on the extension of the ERSEM
approach to the system state variables of the
plankton ecosystem. As shown in Tab. 1.1,
each variable is mathematically expressed by a
multi-dimensional array that contains the concentrations of reference chemical constituents. We
use a superscript indicating the CFF for a specific
living functional group and a subscript for the basic constituent. For instance, diatoms are LFG of
producers and comprise 6 living CFFs written as


(1)
(1) (1) (1) (1) (1) (1)
.
Pi ≡ Pc , Pn , Pp , Ps , Pf , Pl
The standard model resolves 4 different groups
for phytoplankton P( j) , j = 1, 2, 3, 4 (diatoms, autotrophic nanoflagellates, picophytoplankton and
large phytoplankton), 4 for zooplankton Z ( j) , j =
3, 4, 5, 6 (carnivorous and omnivorous mesozooplankton, microzooplankton and heterotrophic
nanoflagellates), 1 for bacteria, 7 inorganic variables for nutrients and gases (phosphate, nitrate,
ammonium, silicate, reduction equivalents, oxygen, carbon dioxide) and 10 organic non-living
components for dissolved and particulate detritus
(cf.. Tab. 1.1 and Fig. 1.1). The state variable nitrate is assumed here to be the sum of both nitrate
and nitrite. Reduction equivalents represent all
the reduced ions produced under anaerobic conditions. This variable was originally used only
in the benthic nutrient regeneration module of
ERSEM (Ruardij and Van Raaphorst, 1995) but
was extended to the water column in Vichi et al.
(2004). In this approach, all the nutrient-carbon
and chlorophyll-carbon ratios are allowed to vary

within their given ranges and each component has
a distinct biological time rate of change.
Following Vichi et al. (2007b) the dynamical
equations are written in two forms: 1) rates of
change form; and 2) explicit functional form. In
“rates of change form”, the biogeochemical reaction term a generic state variable C is written as:
dC
dt

=
bio

∑ ∑

i=1,n j=1,m

dC
dt

ej

,

(1.0.1)

Vi

where the right hand side contains the terms representing significant processes for each living or
non-living component. The superscripts e j are
the abbreviations indicating the process which determines the variation. In Table 1.2 we report
the acronyms of the processes used in the superscripts. The subscripts Vi indicate the state variable involved in the process. If V = C, we refer to
intra-group interactions such as cannibalism.
When a term is present as a source in one equation and as a sink in another, we refer to it following this equivalent notation:
dC
dt

e

=−
V

dV
dt

e

.

(1.0.2)

C

In “functional process form”, the formulation
of the dynamic dependencies on other variables
is made explicit, i.e.: all the rates of change in
eq. (1.0.1) are given in the complete functional
parameterization. Although this is the more complete mathematical form, it is more difficult to
read and interpret at a glance, especially when trying to distinguish which processes affect the dynamics of which variable . All equations in this
Part will be presented both in rate of change and
in functional process forms.

15

1. The formalism of the BFM equations

Figure 1.1.: Scheme of the state variables and pelagic interactions of the reference STANDALONE
model. Living components are indicated with bold-line square boxes, non-living organic
components with thin-line square boxes and inorganic components with rounded boxes
(modified after Blackford and Radford (1995) and Vichi et al. (2007b)).

16

Variable

Code

Type

Const.

Units

Description

N (1)

N1p

IO

P

mmol P m−3
m−3

Phosphate

N (3)

N3n

IO

N

mmol N

N (4)

N4n

IO

N

mmol N m−3

N (5)

N5s

IO

Si

mmol Si m−3

Silicate

N (6)

N6r

IO

R

mmol S m−3

Reduction equivalents, HS−

N (7)
O(2)

N7f
O2o

IO
IO

Fe
O

Dissolved Iron
Dissolved Oxygen

O(3)

O3c

IO

C

O(5)

O3h

IO

-

µ mol Fe m−3
mmol O2 m−3
mg C m−3
mmol Eq m−3

(1)
Pi

P1[cnpslf]

LO

m−3 ,

Nitrate
Ammonium

C N P Si

mg C

Chl Fe

N-P-Si m−3 , mg

mmol

Dissolved Inorganic Carbon
Total Alkalinity
Diatoms

Chl-a m−3 , µ mol
Fe m−3
(2)

Pi

P2[cnplf]

LO

C N P Chl

mg C m−3 , mmol

Fe

N-P m−3 , mg Chl-a

NanoFlagellates

m−3 , µ mol Fe m−3
(3)
Pi

P3[cnplf]

LO

C N P Chl

mg C m−3 , mmol

Fe

N-P m−3 , mg Chl-a

Picophytoplankton

m−3 , µ mol Fe m−3
(4)
Pi

P4[cnplf]

LO

C N P Chl

mg C m−3 , mmol

Fe

N-P m−3 , mg Chl-a

Large phytoplankton

m−3 , µ mol Fe m−3
Bi

B1[cnp]

LO

CNP

mg C m−3 , mmol
N-P

(3)

Zi

Z3[cnp]

LO

CNP

mg C m−3 , mmol
N-P

(4)

Zi

Z4[cnp]

LO

CNP

Z5[cnp]

LO

CNP

(1)

Ri

Omnivorous Mesozooplankton

m−3

mg C m−3 , mmol
N-P

(6)

Zi

Carnivorous Mesozooplankton

m−3

mg C m−3 , mmol
N-P

(5)

Zi

Pelagic Bacteria

m−3

Microzooplankton

m−3

Z6[cnp]

LO

CNP

mg C m−3 , mmol
N-P m−3

Heterotrophic Flagellates

R1[cnpf]

NO

C N P Fe

mg C m−3 , mmol

Labile Dissolved Organic Matter

N-P m−3 , µ mol Fe
m−3
(2)
Rc
(3)
Ri
(6)
Ri

R2c

NO

C

mg C m−3
m−3

R3c

NO

C

mg C

R6[cnpsf]

NO

C N P Si Fe

mg C m−3 , mmol

Semi-labile Dissolved Organic Carbon
Semi-refractory Dissolved Organic Carbon
Particulate Organic Detritus

N-P-Si m−3 , µ mol
Fe m−3

Table 1.1.: List of the reference state variables for the pelagic model. Type legend: IO = Inorganic; LO
= Living organic; NO = Non-living 
organic. The subscript i indicates
 the basic components
(1)
(1) (1) (1) (1) (1) (1)
(if any) of the variable, e.g. Pi ≡ Pc , Pn , Pp , Ps , Pl , Pf .

17

1. The formalism of the BFM equations
Abbreviation
gpp
rsp
prd
rel
exu
lys
syn
nit/denit
scv
rmn

Process
Gross primary production
Respiration
Predation
Biological release: Egestion, Excretion
Exudation
Lysis
Biochemical synthesis
Nitrification, denitrification
Scavenging
Biochemical remineralization

Table 1.2.: List of all the abbreviations used to indicate the physiological and ecological processes in
the equations

18

2. The pelagic plankton model
2.1. The environmental
parameters affecting
biological rates
This section describes the dependencies of the
pelagic biogeochemical processes from the physical environment. The BFM has defined a set of
environmental variables that are provided by external data or by a physical model (Tab.2.1). The
local coupling between physics and biogeochemistry is realized explicitly through surface irradiance and temperature. Temperature regulates all
physiological processes in the model and its effect, denoted by f T , is parameterized in this nondimensional form
T −10

f T = Q1010

(2.1.1)

where the Q10 coefficient is different for each
functional process considered.
Light is fundamental for primary producers
and the energy source for photosynthesis is the
downwelling amount of the incident solar radiation at the sea surface. Only a portion of the
downwelling irradiance is used for growth, the
Photosynthetic Available Radiation (PAR) EPAR
(the notation of Sakshaug et al. (1997) is used
here). The BFM only computes the propagation
of PAR as a function of the attenuation coefficient of sea water, which is determined by the
concentration of suspended and dissolved components. The BFM considers two cases controlled by the parameters described in Tab.5.3: 1)
a broadband mean attenuation coefficient and 2) a
3-band attenuation coefficient for chlorophyll as
proposed by Morel (1988) and tabulated according to Lengaigne et al. (2007).
In the default case we assume that PAR is propagated according to the Lambert-Beer formulation with broadband, depth-dependent extinction
coefficients
EPAR (z) = εPAR QS eλw z+

R0
z

λbio (z′ )dz′

where εPAR is the coefficient determining the fraction of PAR in QS . Light propagation takes into
account the extinction due to suspended particles,
λbio , and λw as the background extinction of water. The biological extinction is written as
4

λbio =

∑c

j=1

P( j)

( j)

Pl + c

R(6)

(6)

Rc

(2.1.3)

where the extinctions due to the concentration of
phytoplankton chlorophyll and particulate detritus are considered (see Sec. 2.2 and 2.5). The c
constants are the specific absorption coefficients
of each suspended substance (Tab. 2.2 and 2.8),
which may also include extinction by suspended
sediments (not currently resolved by the BFM but
can be included from an external model).
With the 4-band wavelength parameterization
the model considers differential attenuation of
light divided in red, green and blue coefficients
(RGB) that are functions of the chlorophyll concentration. The visible portion is divided equally
in the 3 bands such as
 R0

R0
R0
′
′
′
′
ε
′
′
EPAR (z) = PAR QS e z λR (Chl,z )dz + e z λG (Chl,z )dz + e z λB (Chl,z )dz
3
(2.1.4)
where the λ coefficients are obtained with a lookup table dependng on the local total chlorophyll
( j)
concentration Chl = ∑4j=1 Pl . In this parameterization both the background attenuation of water
and the specific attenuation of chlorophyll are not
used.
The short-wave surface irradiance flux QS
is obtained generally from data or from an
atmospheric radiative transfer model. Three
types of forcings are allowed depending
on the value of LightPeriodFlag in
Pelagic_Environment.nml:

(2.1.2)

19

2. The pelagic plankton model
1. Instantaneous
(LightPeriodFlag=1)

irradiance

2. Daily average (LightPeriodFlag=2)
3. Daylight average (with explicit photoperiod,
LightPeriodFlag=3)
k
The discretized irradiance EPAR
is located
at the top of each layer k, and three different choices are available according to
the parameter LightLocationFlag in
Pelagic_Ecology.nml:
k
1. light at the upper interface EPAR = EPAR
(LightLocationFlag=1),

2. light
in
the
middle
of
the
cell
as
in
Lazzari et al.
(2012)
(LightLocationFlag=2),


k
k ∆z
k
(2.1.5)
EPAR = EPAR exp −λ
2
3. integrated
light
over
the
level
depth
as
in
Vichi et al.
(2007b)
(LightLocationFlag=3),

(2)

• autotrophic nanoflagellates (Pi ), ESD =
2-20µ , motile unicellular eukaryotes comprising smaller dinoflagellates and other autotrophic microplanktonic flagellates eaten
by heterotrophic nanoflagellates, micro- and
mesozooplankton;
(3)

• picophytoplankton (Pi ), ESD = 0.2-2µ ,
prokaryotic organism generally indicated
as non-diazotrophic autotrophic bacteria
such as Prochlorococcus and Synechococcus, but also mixed eukaryotic species
(Worden et al., 2004), mostly preyed by heterotrophic nanoflagellates;
(4)

• large, slow-growing phytoplankton (Pi ),
ESD > 100µ , that represents a wide group
of phytoplankton species, also comprising
larger species belonging to the previous
groups (for instance dinoflagellates) but also
those that during some period of the year develop a form of (chemo)defense to predator attack. This group generally has low
growth rates and small food matrix values
with respect to micro- and mesozooplankton
groups.

The processes parameterized in the biological



k
EPAR
k k
source term are gross primary production (gpp),
∆z
λ
1
−
exp
−
λ k ∆zk
(2.1.6) respiration (rsp), exudation (exu), cell lysis (lys),
nutrient uptake (upt), predation (prd) and bioWhen needed for physiological computations, the chemical synthesis (syn). All the phytoplankton
value is converted from W m−2 to the unit of groups share the same form of primitive equaµ E m−2 s−1 with the constant factor 1/0.215 tions, but are differentiated in terms of the values
(Reinart et al., 1998).
of the physiological parameters.
Phytoplankton in the model is composed of
several constituents and the related dynamical
2.2. Phytoplankton
equations (Box 2.1): some of them are mandaPrimary producers are divided in four principal tory (C, N, P), that is they are required to compute
sub-groups coarsely representing the functional physiological terms, and some others are required
spectrum of phytoplankton in marine systems. for specific groups (like Si in diatoms) or they
The operational model definitions of the phyto- are optional such as the Chl and Fe constituents.
The inclusion of Fe is a typical example of the
plankton groups are:
model flexibility as it allows to have an additional
(1)
• diatoms (Pi ), ESD = 20-200µ , unicellular multiple-nutrient limitation term without affecteukaryotes enclosed by a silica frustule eaten ing the other parameterizations (see Sec. 2.2.8).
When Chl constituent is not used (see Sec.
by micro- and mesozooplankton;
2.2.7), light acclimation in phytoplankton can be
EPAR =

20

2.2. Phytoplankton

Code

Unit

Description

T

ETW

oC

Sea water salinity

S

ESW

(-)

Sea water temperature

ρ
EPAR
SS
Watm
Fice
Φ

ERHO

Kg m−3

Sea water density

EIR

µ E m−2 s−1
g m−3
m s−1
(-)
hours

Photosynthetically available radiation

Symbol

ESS
EWIND
EICE
SUNQ

Suspended sediment load
Wind speed
Sea ice fraction
Daylight length

Table 2.1.: Mathematical and code symbols, units and description of the environmental state variables
defined in the model.

Equation Box 2.1 Phytoplankton equations
dPc
dt
dPn
dt
dPp
dt
dPl
dt

dPs
dt

=
bio

=
bio

dPc
dt

∑

i=3,4

=
bio

=
bio

=
bio

gpp

−
O(3)

upt

dPl
dt

syn

dPs
dt

upt

(2)

Rc

−
N (i)

j=1,6

dPp
dt

lys

Pl
dPc
∑
Pc j dt

prd

−
N (1)

∑

j=1,6

−

−
N (5)

dPs
dt

(6)

=
bio

dPf
dt

upt

−
N (7)

dPf
dt

−

j=1,6

dPc
dt

lys
( j)

Pp
dPc
∑
Pc k=4,5,6 dt

−

Rc

Pn
dPc
∑
Pc k=4,5,6 dt

∑

k=4,5,6

dPc
dt

prd
(k)

(2.2.1a)

Zc

prd
(k)

(2.2.1b)

Zc

prd
(k)

(2.2.1c)

Zc

(2.2.1d)

( j)

Rs

Ps
dPc
∑
Pc k=4,5,6 dt

dPs
dt

lys

−

Rf

if Pf ∃, otherwise

( j)

−

∑

Zc

−

(6)

O(3)

Rn

Rp

lys

if Ps ∃, otherwise
dPf
dt

(i)

−

lys

dPn
dt

∑

rsp

dPc
dt

−

upt

dPn
dt

dPp
dt

exu

dPc
dt

dPf
dt

prd
(k)

(2.2.1e)

Zc

=0
bio

Pf
Pc

∑
j

dPc
dt

=0

prd
( j)

(2.2.1f)

Zc

(2.2.1g)

bio

21

2. The pelagic plankton model
described by means of the optimal light property
dPc gpp
as proposed by Ebenhöh et al. (1997). The no= fPT fPE fPPP rP0 Pc .
(2.2.6)
dt O(3)
tation remains mostly unchanged but in this case
the variable P̂l indicates the optimal light level to
Note that it is possible to further control temwhich phytoplankton is acclimated and the alter- perature limitation with a cut-off value cTP applied
native equation to (2.2.1d) is given in Sec. 2.2.7. to the temperature regulating factor as done in
(Vichi et al., 2007a) to control the growth of picophytoplankton at high latitudes:
2.2.1. Photosynthesis and carbon
dynamics


f T = max 0, f T − cTP .
(2.2.7)
Photosynthesis is primarily controlled by light
by means of the non-dimensional light regulation Respiration is defined as the sum of the basal resfactor proposed by Jassby and Platt (1976)
piration, which is independent of the production
rate, and the activity respiration:


EPAR
E
(2.2.2)
fP = 1 − exp −
dPc gpp
dPc rsp
EK
T
= bP fP Pc + γP
.
(2.2.8)
dt O(3)
dt O(3)
where E is the optimal irradiance, defined as
K

Ek = Pm∗ /α ∗
Pc
Pl

Pm∗ =

fPT fPPP rP0

α∗ =

0
fPT fPPP αchl

(2.2.3)

The maximum chl-specific photosynthetic rate Pm∗
is controlled by the non-storable nutrients that
control primary production fPPP (Sec. 2.2.2). It
is assumed that both α ∗ and Pm∗ are controlled
by the same regulating factors (Jassby and Platt,
1976; Behrenfeld et al., 2004). When instantaneous light is used, the equation is written in the
following form:


α 0 EPARPl
fPE = 1 − exp − chl 0
rP Pc

(2.2.4)

while in case the length of the photoperiod is considered, this is implicitly done in Pm∗ that becomes.
Pm∗ = fPT fPPP rP0

Pc Φ
Pl 24

(2.2.5)

where Φ is the daylight length in hours
(Lazzari et al., 2012). The final form of gross primary production also depends on the choice of the
daylight availability. In the case of instantaneous
or average daily irradiance (Vichi et al., 2007a), it
is written as

22

Basal respiration is only a function of the carbon biomass, temperature (through the regulating
factor fPT ) and the specific constant rate bP . The
activity respiration is a constant fraction (γP ) of
the total gross primary production.
The term lysis includes all the non-resolved
mortality processes that disrupt the cell membrane, such as mechanical causes, virus and
yeasts. It is assumed that the lysis rate is partitioned between particulate and dissolved detritus
and tends towards the maximum specific rate dP0
with a saturation function of the nutrient stress as:


hPp,n
n,p
dPc lys
lys
0
= εP
d Pc + χ (2.2.9)
dt R(6)
fPp,n + hPp,n P
c

n,p
dPc lys 
= 1 − εP
(2.2.10)
dt R(1)
c


hPp,n
0
lys
d Pc + χ
fPp,n + hPp,n P
The lysis of cells generates both dissolved and
particulate detritus; the structural parts of the cell
are not as easily degradable as the cytoplasm,
therefore the percentage going to DOC is inversely proportional to the internal nutrient content and constrained by the minimum structural
content in the following way:
!
min
n
pmin
n,p
εP = min 1, P , P
.
(2.2.11)
Pp /Pc Pn /Pc

2.2. Phytoplankton
Symbol

Code

Unit

Description

Q10P

p_q10

-

Characteristic Q10 coefficient

cT
P

p_temp

-

Cut-off threshold for temperature regulating factor

r 0P

p_sum

d−1

Maximum specific photosynthetic rate

bP

p_srs

d−1

Basal specific respiration rate

d0P

p_sdmo

d−1

Maximum specific nutrient-stress lysis rate

hPp,n,s

p_thdo

-

Nutrient stress threshold

p_seo

d−1

Extra lysis rate

hP

p_sheo

mg C m−3

Half saturation constant for extra lysis

βP

p_pu_ea

-

Excreted fraction of primary production

γP

p_pu_ra

-

Activity respiration fraction

p_switchDOC

1-3

Switch for the type of DOC excretion. Choice consistent with bacteria parame-

x

dP
x

(1)

terization in Sec. 2.3. 1. All DOC is released as Rc (BACT1); 2. Activity DOC
(2)

(2)

is released as Rc (BACT2); 3. All DOC is released as Rc (BACT3)
p_netgrowth

logical

Switch for parameterization of nutrient limited growth (T or F)

p_limnut

0-2

Switch for parameterization of nutrient co-limitation

a4

p_qun

m3 mg C−1 d−1

Specific affinity constant for N

n

p_lN4

mmol

Half saturation constant for ammonium uptake preference over ammonium

hP

N-NH4 m−3
, nopt
, nmax
nmin
P
P
P

mmolN mgC−1

Minimum, optimal and maximum nitrogen quota

p_qup

m3 mg C−1 d−1

Specific affinity constant for P

p_qplc,

mmolP mgC−1

Minimum, optimal and maximum phosphorus quota

p_switchSi

1,2

Switch for parameterization of silicate limitation (1=external, 2=internal)

p_chPs

mmol Si m−3

Half saturation constant for Si-limitation
Variable half saturation constant for Si-limitation of Contois (like Monod if 0)

p_qnlc,
p_qncPPY,
p_xqn*p_qncPPY

a1
pmin ,
P

popt ,
P

pmax
P

p_qpcPPY,
p_xqp*p_qpcPPY
hsP

ρPs

p_Contois

-

a5

p_qus

m3 mg C−1 d−1

p_qslc

C−1

smin ,sopt
P

P

mmolSi mg

Specific affinity constant for Si
Minimum and optimal Si:C ratio in silicifiers

p_qscPPY
sink

ωP

sink

lP

p_res

m d−1

p_esNI

-

Maximum sedimentation rate
Nutrient stress threshold for sinking

p_switchChl

Choice of the Chl synthesis parameterization

dPchl

p_sdchl

d−1

Chlorophyll turn-over rate

0
αchl

p_alpha_chl

mgC (mg chl)−1

Maximum light utilization coefficient

µ E−1 m2
0
θchl

p_qlcPPY

mg chl mg C−1

Maximum chl:C quotum

cP

p_epsChla

m2 (mg chl)−1

Chl-specific light absorption coefficient

p_EpEk_or

-

Optimal value of EPAR /EK

p_tochl_relt

d−1

Relaxation rate toward the optimal Chl:C value

p_iswLtyp

0-6

Shape of the productivity function (diagnostic Chl, ChlDynamicsFlag=1)

p_chELiPPY

W m−2

Maximum and minimum thresholds for light adaptation

d−1

Relaxation rate toward the optimal light

ε

opt

τ opt
EPmax ,EPmin

p_clELiPPY
vIP

p_ruELiPPY
p_rPIm

m

d−1

Additional background sinking rate

Table 2.2.: Mathematical and code symbols, units and description of the phytoplankton parameters
(namelist Pelagic_Ecology.nml: Phyto_parameters).
23

2. The pelagic plankton model
This equation ensures that the carbon and nutrients in the supportive structures, which are asmin
sumed to have pmin
P and nP nutrient ratios, are
always released as particulate components.
The model also implements and extra lysis rate
in the last term of eq. (2.2.9-2.2.10), which is dependent on the population density. This term is
controlled by an optional specific lysis rate
Pc
Pc
(2.2.12)
χ lys = dPx
Pc + hxP
which is usually included for those phytoplankton
population that are indedible and therefore acts as
an additional mortality term to regulate the population dynamics.

2.2.2. Multiple nutrient limitation
The BFM inherits the treatment of nutrient limitation from the original work by
Baretta-Bekker et al. (1997) which has similarities to the Caperon and Meyer equation, with
some specific extensions. Nutrients are divided
in two main groups, the ones that directly control carbon phosynthesis (as silicate, since duplication can only occur with Si deposition) and the
ones that are decoupled from carbon uptake because of the existence of cellular storage capabilities. Limiting factors for nutrients can be both internal (i.e. based on the internal nutrient quota) or
external (based on dissolved inorganic concentration). The internal limitation is only partly based
on Droop (1973), with the concepts of nutrient
surge and storage by Baretta-Bekker et al. (1997)
and the possibility to apply the parameterization
of net growth described in Vichi et al. (2004) by
means of the parameter p_netgrowth in the
namelist (Tab. 2.2). The regulating factors for
internal limitation are controlled by the optimal
and minimum values for the existence of strucmin
tural components (nopt
P , nP ), and they are implemented for N, P and Fe (Sec. 2.2.8) and also for
Si (Lazzari et al., 2012), although in this case it is
not possible to have luxury uptake because there
is no storage for dissolved silicate in the cell:
fPn

24





Pn /Pc − nmin
P
= min 1, max 0, opt
nP − nmin
P



(2.2.13)

fPp





Pp /Pc − pmin
P
= min 1, max 0, opt
pP − pmin
P

fPf =

Pf /Pc − φPmin
φPopt − φPmin

fˆPs =



(2.2.14)

(2.2.15)

Ps /Pc − smin
P

(2.2.16)

sopt
− smin
P
P




Pp /Pc − pmin
P
fPp = min 1, max 0, opt
pP − pmin
P
(2.2.17)
Note that the ratios between brackets can be larger
than 1 because the nutrient quota are allowed to
reach the maximum values (cf. Sec. 2.2.3) but
the regulating factor has to be limited to 1 because
above the optimal quotum there is no limitation to
physiological processes.
There is a flag (Tab. 2.2) for the choice of
external and internal silicon limitation in case
of diatoms. External dissolved silicate control
should be preferable (Flynn, 2003) and it implements a Monod regulation with variable halfsaturation constant hsP , modified according to
Contois (1959):
fPs (1) =
fPs ( j)

N (5)
(1)

N (5) + (hsP + ρPs Ps )
= 1, j 6= 1

(2.2.18)

The Contois formulation implements a control on
dissolved nutrient uptake that is dependent on the
size of the population, in this case it is applied
to biogenic silica as suggested by Tsiaras et al.
(2008). If the Contois parameter is zero a classical Monod function is obtained, while if ρPs 6= 0
the half-saturation constant is increased and the
resulting silicate control factor is smaller.
Multiple nutrient limitation is different for nutrients that can be stored in the cell and nutrients
that cannot. The threshold combination (Liebiglike) is usually preferred to the multiplicative approach (Flynn, 2003): the BFM allows the three

2.2. Phytoplankton
alternative ways implemented in ERSEM-II to only structural), and a rate based on consideracombine N and P limitation (that can be selected tions of balanced growth and luxury uptake:
differently for each phytoplankton group, Tab.
dPn upt
2.2),
∑ dt ( j) =
N
j=3,4




fPn,p = min fPn , fPp
(2.2.19)
hnP
n
(3)
n (4)
min
aP n
N + aP N
Pc ,
2
hP + N (4)
(2.2.20)
fPn,p =
p
n

 
1/ f + 1/ fP
Pn
q P
max
opt
nP GP + νP nP −
Pc
(2.2.23)
fPn,p =
fPn fPp
(2.2.21)
Pc

A preference for ammonium is parameterized
through a saturation function that modulates the
affinity constant for dissolved nitrate. Balanced
uptake is a function of the net primary production
G p , which is computed as:
!
dPc gpp dPc exu dPc rsp dPc lys
GP = max 0,
−
−
−
dt O(3) dt R(i)
dt O(3) dt R(i)
c
c
(2.2.24)
and


GP
and it is considered in the parameterization of
.
νP = max 0.05,
some processes such as chlorophyll synthesis and
Pc
sinking (Secs. 2.2.6 and 2.2.9).
When the nitrogen uptake rate in eq. (2.2.23)

and a simple multiplicative approach for the effect of nutrients that directly control photosynthesis (eq. 2.2.6), such as f PP = fˆPs as done by
Lazzari et al. (2012) or f PP = fPf fPs when iron dynamics is included (Sec. 2.3) as in Vichi et al.
(2007b). The co-limitation from all nutrients is
always done with a threshold method


fPnut = min fPn,p , fPf , fPs
(2.2.22)

is positive, the partitioning between N (3) and
2.2.3. Nutrient uptake
N (4) uptake is done by multiplying the rate by the
A major problem in modelling nutrient up- fractions:
hnP
take is to deal with unbalanced growth condi(3)
anP hn +N
(4) N
(3)
P
tions, and any realistic application is expected
(2.2.25)
εP =
hnP
(3)
anP N (4) + anP hn +N
to produce transient environmental conditions re(4) N
P
sulting in uncoupled assimilation rates of carbon and nutrients. The Droop (intracellular
anP N (4)
(4)
=
(2.2.26)
ε
P
hnP
quota approach) and Monod (external concentra(3)
anP N (4) + anP hn +N
N
(4)
P
tion approach) equations produce consistent results when applied in balanced growth condi- When it is negative, the whole flux is directed to
tions (Morel, 1987). The BFM approach, largely the DON pool R(1)
n .
derived from Baretta-Bekker et al. (1997), combines both mechanisms with a threshold control. 2.2.3.2. Phosphorus
2.2.3.1. Nitrogen
The uptake of DIN is the sum of the uptake of dissolved nitrate and ammonium and is the minimum
between a diffusion-dependent uptake rate (when
internal nutrient quota are low and nutrients are

Phosphate uptake is simpler than N uptake as one
single species is considered:
dPp
dt

upt
N (1)


(1)
= min aP N (1) Pc , popt
P GP +

 
Pp
max
+νP pP −
Pc
(2.2.27)
Pc

25

2. The pelagic plankton model
2.2.3.3. Silicate
Silicate is not stored in the cell and therefore the
uptake is directly proportional to the net carbon
growth (see eq. 2.2.24):
dPs
dt

upt
N (5)

= sopt
P GP

(2.2.28)

2.2.4. Nutrient loss associated with
lysis

of DOC (Tab. 1.1) and phytoplankton excretion
must be consistent with the bacteria parameterization (Sec. 2.3). Carbohydrates are excreted
when phytoplankton cannot equilibrate the fixed
C with sufficient nutrients to maintain the minimum quotum needed for the synthesis of new
biomass. The choice of the parametrization for
DOC exudation is controlled by the flag parameters p_netgrowth and p_switchDOC in the
phytoplankton namelist (Tab. 2.2), which can
be set differently for each sub-group. There
is a runtime check that controls the consistency between the two parameters. When the
flag is false, the ERSEM-II parameterization
(Baretta-Bekker et al., 1997) as implemented in
Vichi et al. (2007b) is used,

Nutrients in the cell are not evenly distributed between cytoplasm and structural parts. When the
cell wall is disrupted, part of the nutrients are released as dissolved organic matter (the nutrients
in the cytoplasm) and part as particulate organic
matter (the nutrients in the structural parts, such
as the cell wall).

 dPc gpp
dPc exu 
= βP + (1 − βP ) 1 − fPn,p
The redistribution of the cellular elements C,
dt R(1)
dt O(3)
c
N and P over the excretion variables reflects the
(2.2.32)
preferential remineralization of P and to a lesser
where exudation (and thus carbon biomass
extent of N, with respect to C. This parameter valloss) increases when phytoplankton has low nuues are based on the ERSEM-II parameterizations
trient:carbon ratios.
The nutrient-stress ex(Baretta-Bekker et al., 1997).
cretion can also be partly directed to semi(2)
refractory DOC (Rc ) by setting the parameter
lys
p,n
dPp
h
p_switchDOC=3:
= εPi p,n P p,n dP0 Pi i = n, p (2.2.29)
dt R(6)
fP + hP
i
dPc exu
dPc gpp
=
(2.2.33)
β
P
p,n
dt R(1)
dt O(3)


dPp lys
h
c
P
= 1 − εPi
d 0 Pi i = n, p

 dPc gpp
dPc exu
dt R(1)
fPp,n + hPp,n P
i
= (1 − βP ) 1 − fPn,p
(2.2.34)
dt R(2)
dt O(3)
(2.2.30)
c
In the case of Si, the release is always in particulate form and thus is a constant fraction of the However, the ERSEM-II parameterization not
only leads to the desired extra release of Cparticulate carbon lysis
enriched DOM but also to a decrease of the population biomass. It has been shown by Vichi et al.
(1) lys
(1) lys
(2004) that the use of a different parameterizadPc
dPs
= sopt
(2.2.31)
P
tion results in the maintenance of higher standdt
dt
(6)
(6)
Rs
Rc
ing stocks in the post bloom period with a consequently enhanced consumption of the non2.2.5. Exudation of carbohydrates
limiting nutrients. This parameterization is also
In the case of intra-cellular nutrient shortage, used by Lazzari et al. (2012) to allow the connot all photosynthesized carbon can be assimi- sumption of N species in a P-depleted environlated into biomass, and the non-assimilated part ment.
When the p_netgrowth flag is true, total exis released in the form of dissolved carbohy(2)
drates. The model considers three different types udation is directed to the state variable Rc (see

26

2.2. Phytoplankton
also Sec. 2.3 for a detailed explanation of DOM thesis:
in the BFM), and it is written as the sum of an activity exudation linked to the gross primary production and a balance term

ρchl =

gpp
dPc
dt O(3)
0
θchl
∗
α EPAR Pl

(2.2.38)

0
and multiplying by a maximum chl:C ratio θchl
+ GP − Gbal
(2.2.35) which is different for each phytoplankton funcP
(2)
(3)
Rc
O
tional group (Tab. 2.2). The major difference in
eq. (2.2.37b) proposed by Lazzari et al. (2012) is
where
the presence of the combined nutrient regulating
factor, which implies that chl synthesis is reduced
upt
in regions limited by phosphorus like for instance
1
dPn
Gbal
,
P = max 0, min GP , min ∑
the Mediterranean.
nP j=3,4 dt N ( j)
Following the notation shown in Sec. 2.2.3,
(2.2.36)
Geider’s
original formulation is rewritten after

1 dPp upt
some algebra as:
dt N (1)
pmin
P
fPE rP0 Pc
0
ρchl = θchl
(2.2.39)
0 E
αchl
PAR Pl
2.2.6. Chlorophyll synthesis and

dPc
dt

exu

dPc
= βP
dt

gpp

photoacclimation
The chlorophyll equation in (2.2.1d) is composed
of two terms. The first one is net chlorophyll synthesis, which is mostly derived from Geider et al.
(1996, 1997) with some adaptations, and the second one represents the losses due to grazing.
Net chlorophyll synthesis is namely a function
of acclimation to light conditions, nutrient availability and turnover rate. The former process is
taken into account by Geider’s parameterization,
while the latter is generally parameterized with
different formulations, for instance by assuming a
dependence on gross carbon uptake (Geider et al.,
1997; Blackford et al., 2004) and/or on nitrogen assimilation (Geider et al., 1998; Flynn et al.,
2001). To integrate the variety of processes
into our formulations, we consider three different parameterizations presented in 2.2 that have
been used in different applications of the BFM
(Vichi and Masina (2009), Lazzari et al. (2012),
Clementi et al., in preparation). The synthesis
part is rather similar as it is controlled by the dynamical chl:C ratio ρchl proposed by Geider et al.
(1997), which regulates the amount of chlorophyll in the cell according to a non-dimensional
ratio between the realized photosynthetic rate in
eq. (2.2.6) and the maximum potential photosyn-

The ratio is down-regulated when the rate of light
absorption (governed by the quantum efficiency
and the amount of light-harvesting pigments) exceeds the rate of utilization of photons for carbon fixation, as explained in detail in Geider et al.
(1996). Also here it is assumed that both α ∗ and
Pm∗ are controlled by the same regulating factors
(Jassby and Platt, 1976; Behrenfeld et al., 2004).
The formulation of the lysis term is still unknown and this is reflected in the diversity
of parameterizations that are available in the
model and controlled by the namelist parameter
p_switchChl (Tab. 2.2). The user should consider the different assumptions and use them accordingly. A different parameterization can be
used for each phytoplankton group to allow their
testing. In eq. (2.2.37a) it is assumed that chl lysis
is simply linked to the carbon term, while in eq.
(2.2.37b) there is a relaxation to the optimal chl:C
ratio with an additional lysis term linked to background mortality. The parameterization proposed
in eq. (2.2.37c) is slightly more elaborated as it
considers an optimal cell acclimation to light that
corresponds to an optimal value EPAR/EK = ε opt
of the exponential in eq. (2.2.4).
The theoretical chlorophyll concentration corresponding to optimal light acclimation in eq.
(2.2.37c) is computed as follows:

27

2. The pelagic plankton model
Equation Box 2.2 Different parameterizations for the chlorophyll synthesis.
dPl
dt
dPl
dt

syn

dPl
dt

syn

dPl
dt

syn

syn

Pl dPc
= ρchl GP −
Pc dt

lys

(2.2.37a)



fPp,n ρchl GP − max 0, dP◦ (1 − fPp,n ) Pl −


0
min (0, GP ) ∗ max 0, Pl − θchl
Pc
!


dPc lys dPc rsp
= ρchl GP − θchl
+
− max 0, Pl − Plopt τchl
dt
dt
=

rP0 Pc
0
Plopt αchl
Plopt

=

=

ρchl GP − dPchl Pl

p,n,s

hP

EPAR
ε opt

= ε opt

rP0 Pc
0 E
αchl
PAR

(2.2.40)

and the chl content in (2.2.37c) is relaxed to this
optimal value with a specific time scale parameter
τchl .
The parameterization presented in eq.
(2.2.37d) has instead an explicit term for
the turn-over rate of chl, dPchl . The losses of
chlorophyll are not considered as mass losses in
the model because we have currently not implemented a chl component in detritus and dissolved
organic matter. The same consideration applies to
the ingested chl fraction in zooplankton. All these
terms are presently collected into a generic sink
term used for mass conservation purposes, which
can be easily split into its major components once
it is deemed necessary to follow the degradation
products of chl (phaeopigments).

2.2.7. Light limitation and
photoacclimation based on
optimal light
This parameterization is alternative to the one
above and is consistent with the one originally implemented in ERSEM-II (Ebenhöh et al., 1997).
It assumes that phytoplankton is acclimated to the
prevailing light in a few days and that the C:chl
ratio is constant.
28

−

fPnut




1
dPc
Pl
−
Pc 1 + EPAR dt

(2.2.37b)

(2.2.37c)

rsp

(2.2.37d)

The state variable Iopt in the original ERSEMII formulation represented the light intensity at
which production saturated for the whole phytoplankton community. We consider here an optimal light for each phytoplankton group indicated
( j)
with P̂l , which has the same meaning as the
light saturation parameter Ek used in the alternative formulation of chlorophyll dynamics (Sec.
2.2.3).
A note of caution on the application of this
parameterization: the model was formulated to
simulate the daily production, therefore it is suggested to apply it with the daylight-averaged irradiance, possibly modulated by the duration of the
daylight period (see Sec. 2.1). Several forms of
the P-E curve are available:


EPAR
prod = fPT fPs r0P min 1,
(ramp)
P̂l
prod = fPT fPs r0P

EPAR 1− EPAR
e P̂l
P̂l

(Steele)

The light regulating factor is computed by integrating the P-E curve over the depth of the considered layer, therefore the irradiance at the top of
the layer is used
1
fP = T s
fP fP r0P D
E

Z 0

−D

prod



EPAR
P̂l



dz (2.2.41)

2.3. Bacterioplankton
The time rate of change of P̂l is computed with
a relaxation term
 opt

d P̂l
(2.2.42)
= νPI EP − P̂l
dt bio

from dissolved pools is computed as for N and P
(Sec. 2.2.3) by taking the minimum of two rates, a
linear function of the ambient concentration simulating the membrane through-flow at low external
Fe concentration, and the balancing flux accordopt
where EP is the reference light saturation param- ing to the carbon assimilation:
eter to which the phytoplankton adapt with frequency νPI (generally 4 days for a full acclima
∂ Pf upt
tion).
=min a7P N (7) Pc , φPopt GP + (2.2.43)
∂ t N (7)
The irradiance level to which phytoplankton

 
may adapt is
Pf
T 0
max
+ fP rP φP −
Pc
Pc



opt
EP = min EPmax , max EPmin , EPAR
It is assumed that the only physiological iron loss
which is a constrained ramp function over the from phytoplankton is linked to cell disruption,
computed according to carbon lysis and assuming
range of saturating irradiance.
that particulate material has the minimum structural Fe:C ratio:
2.2.8. Iron in phytoplankton
The iron (Fe) constituent in phytoplankton shown
in eq. (2.2.1f) is not activated by default in
the model. It is available as a compilation
key INCLUDE_PELFE and implements the iron
dynamics proposed by Vichi et al. (2007b) described by the parameters of Tab. 2.3. The iron
cycle involves phytoplankton, particulate and dissolved component (Sec. ), while the zooplankton
and bacteria fraction is neglected.
The equation (2.2.1f) for iron in phytoplankton Pf contains a term for the uptake of Fe, a
loss term related to turnover/cell lysis and a predation term. Similarly to N and P content, intracellular Fe:C quota are allowed to vary between
a maximum and a minimum thresholds (φPmax and
φPmin , Tab. 2.3), and the realized quotum is used
to derive a non-dimensional regulating factor as
in eq. (2.2.19). The allowed minimum ratio φPmin
represents the evolutive adaptation of each functional group at the prevailing iron concentrations,
and the optimal value φPopt indicates the cellular
requirement for optimal growth. This regulating
factor modulates the actual photosynthetic rate in
eq. (2.2.6), since there is a clear decrease in the
activity of PSUs due to insufficient cellular Fe
(Sunda and Huntsman, 1997).
The regulating factor inhibits carbon fixation,
but iron can still be taken up in the cell, progressively increasing the internal quotum. Iron uptake

∂ Pf
∂t

lys
(6)

Rf

= φPmin

∂ Pc
∂t

lys
(6)

.

(2.2.44)

Rc

2.2.9. Phytoplankton sinking
The sinking of biogenic material is a fundamental
process for the simulation of carbon sequestration
in the interior of the ocean. However, the estimation of the sinking velocity wB is still parameterized in a very simplified way in the model. Any
phytoplankton group is allowed to have a sinking
velocity using the original ERSEM formulation
(Varela et al., 1995). This is generally valid only
for diatoms, which reach their maximum velocsink
ity ω as a function of the total nutrient stress
(2.2.22) as follows:


sink
(2.2.45)
w (1) = ω max 0, l sink − fPnut
P

where l sink is the nutrient regulating factor value
below which the mechanism is effective.

2.3. Bacterioplankton
Bacterioplankton is included in the standard BFM
with one single state variable representing a wide

29

2. The pelagic plankton model
Symbol

Code

Unit

Description

a7

p_quf

m−3 mg C−1

Specific affinity constant for Fe in phytoplank-

d−1

ton

µ mol
Fe mgC−1

Minimum, optimal and maximum iron quota in

φPmin , φPopt , φPmax

p_qflc,
p_qfcPPY,

phytoplankton

p_xqf*p_qfcPPY
Λ1rmn
f

p_sR1N7

d−1

Specific remineralization rate of dissolved biogenic iron

Λ6rmn
f

p_sR6N7

d−1

Specific remineralization rate of particulate biogenic iron

Q10 f

p_q10R6N7

(-)

Characteristic Q10 coefficient for Fe remineralization

ϕ scv
f
Λscv
f
Λdep
f

m−3

p_N7fsol

µ mol Fe

p_scavN7f

d−1

Specific scavenging rate of dissolved iron

p_qflc

(-)

Specific dissolution fraction of dust iron

Solubility concentration

Table 2.3.: Mathematical and code symbols, units and description of the iron cycle parameters
(namelist
Pelagic_Ecology.nml:Phyto_parameters_iron
Pelagic_environment.nml:PelChem_parameters_iron).
group of aerobic and anaerobic bacteria. The
user can eventually expand the number of bacteria group using the modular facilities described in
Part II of this document and by providing new dynamical equations in the code. The main source of
carbon for bacterioplankton is the organic matter
pool that is composed of particulate detritus (variables R(6) ) and dissolved organic matter (DOM,
variables R(1) and R(2) ).
In addition to the original ERSEM parameterization (Baretta-Bekker et al., 1995) further expanded in Vichi et al. (2007b), the code contains also two alternative parameterizations proposed by Vichi et al. (2004) and Polimene et al.
(2006). We will refer to BACT1 when indicating
the standard most simple parameterization and to
BACT2 and BACT3 for theVichi et al. (2004) and
Polimene et al. (2006) parameterizations, respectively. BACT1 set of equations is the default
choice but all parameterizations are deemed equal
and can be activated at run time.
The different parameterizations consider that
dissolved organic matter (DOM) available to
pelagic bacteria may have different degrees of lability/refractivity (see Tab. 1.1). The lability/refractivity characteristics of DOM are dependent

30

on two factors: the C:N and C:P ratios of bacteria and the structure of organic molecules constituting the DOM matrix. DOM is assumed to be
partitioned into three broad and distinct state variables, each of them corresponding to different degrees of lability/refractivity and having different
production pathways. The parameterizations described below consider only one, two or all of the
3 classes. It is important to remark that the DOC
fraction considered by the model covers only the
labile and the semi-labile DOC part. Therefore, in
the following we define, as “semi-refractory” the
DOC fraction more slowly remineralized by bacteria only to avoid confusion with the truly refractory DOC (with turnover time of 100-1000 years)
that is not considered at all in the model.
The most labile fraction of the total DOM
(1)
pool Ri , i = c, n, p is produced by phytoplankton, zooplankton and bacteria via the lysis, mortality and sloppy feeding (mesozooplankton only)
processes. This DOM variable is characterized by C:N and C:P ratios identical to those
of the producing functional groups. The characteristic turnover time-scale is assumed to be 1
(2)
day. The semi-labile DOM fraction Rc is produced through excretion by phytoplankton and

2.3. Bacterioplankton
bacteria in order to achieve/maintain their internal “optimal” stoichiometry (this part is considered only in BACT2 and BACT3 parameterizations). The production process of semi-labile
DOM can be thought as release of excess carbon and, therefore, negligible N and P pools are
assumed. The characteristic turnover time-scale
is assumed to be 10 days. DOM released by
(3)
bacteria as capsular material Rc is also introduced in BACT2 and BACT3, and it represents
the semi-refractory fraction of the DOM, but it is
only considered as a source of carbon in BACT3
(Polimene et al., 2006). This component of the
DOM pool is also assumed (as for the semilabile fraction) to be DOC only and to be formed
by high molecular weight substances such as
polysaccharid fibrils (Heissenberger et al., 1996),
which are quite resistant to enzymatic attack
(Stoderegger and Herndl, 1998). Therefore, the
characteristic turnover time-scale is assumed to
be longer (100 days), as this material is only
degradable by bacteria at time-scales of 2-3 orders of magnitude longer with respect to the labile
DOC (Stoderegger and Herndl, 1998).
Bacteria have 3 free-varying constituents C, N,
P, and the fundamental equations shown in Eq.
Box 2.3 are common to all parameterizations. In
all versions, pelagic bacteria behave as remineralizers or phytoplankton competitors depending on
their internal nutrient quota, taking up inorganic
nutrients directly from the water. Before providing a detailed description of the three versions the
common processes are described.

metabolism under aerobic conditions and anaerobic metabolism.
The non dimensional regulating factor based on
the nutritional content of bacterial cells is given
by:


B p /Bc Bn /Bc
;
fBn,p = min 1, opt , opt
p
n

(2.3.3)

opt
nopt
B and pB are the “optimal” N : C and P : C
bacterial internal quota (Goldman and McCarthy,
1978).

2.3.2. Respiration
The respiration sink term is composed of basal
and activity respiration as:
dBc
dt

rsp
O(3)

= bB fBT Bc +


a

o

γB + γB (1 −

(2.3.4)
fBo )



dBc
∑ dt
j=1,2,6

upt

( j)

Rc

The basal respiration is parameterized as for
phytoplankton with a constant specific respiration
rate bB and the regulating factor for temperature
given in eq. (2.1.1). γBa and γBO are the fraction of
production that is used for activity respiration under oxic and low oxygen conditions respectively.
Since the BFM pelagic bacteria parameterization
encompasses both aerobic and anaerobic bacterial activities, we consider here the differences in
the energetics of the metabolic pathways in re2.3.1. Regulating factors
lation to oxygen availability. Anaerobic bacteria
have a lower efficiency because they need to conThe oxygen non dimensional regulating factor fBo
sume (respire) more carbon in order to produce
(eq. 2.3.2) is parameterized with a type III control
the same amount of energy.
formulation as:

3
O(2)
2.3.3. Mortality
o
(2.3.2)
fB =

3

3
(2)
o
O
+ hB
The mortality (lysis) process is given by
where the dissolved oxygen concentration O(2) is
considered, and hoB is the oxygen concentration at
which metabolic functionalities are halved. This
steep sigmoid has been proposed by Vichi et al.
(2004) to efficiently switch between bacteria

dBc
dt

lys
(1)

Ri



= d0B fBT + dBd Bc Bc

i = c, n, p
(2.3.5)

31

2. The pelagic plankton model

Symbol

Code

Units

Description

Parameterization

Q10B

p_q10

(-)

Characteristic Q10 coefficient

All

hoB

p_chdo

mmolO2 m−3

Half saturation value for oxygen limitation

All

r0B

p_sum

d−1

Potential specific growth rate

All

bB
γBa

p_srs
p_pu_ra

d−1
(-)

Basal specific respiration rate
Activity respiration fraction

All
All

γBo

p_pu_ra_o

(-)

Additional respiration fraction under anoxic

All

conditions
d0B

p_sd

d−1

Specific mortality rate

All

dd

p_sd2

d−1 mgC−1 m3

Density-dependent specific mortality rate

All

p_suR2

d−1

Specific potential R(2) uptake (semi-labile)

BACT2,BACT3

p_suR3

d−1

Specific potential

p_suR6

d−1

Specific potential R(6) uptake (particulate)

p_suhR1

d−1

R(1)

B

νBR

(2)

(3)
νBR
(6)
νBR
(1)
νBR

R(3)

uptake (semi-refractory)

Specific quality-dependent potential

BACT3
All

up-

All

take (nutrient-rich labile)
(1)

ν0RB

p_sulR1

d−1

Specific quality-independent potential R(1) up-

a1B

p_qup

m−3 mg C−1

Specific affinity constant for P

BACT2

Specific affinity constant for N

BACT2

d−1

Relaxation time scales for N and P uptake or

All
BACT3

BACT1

take (nutrient-poor labile)
d−1

a4B

p_qun

m−3 mg C−1
d−1

νBn , νBp

p_ruen, p_ruep

νBc

p_rec

d−1

remineralization
Relaxation time scale for semi-labile carbon re-

nopt
, popt
B
B

p_qncPBA

mmolN mgC−1 ,

Optimal nutrient quota

All

Minimum nutrient quota

BACT2

Half saturation for nutrient uptake

All

Fractional excretion of R(3) (semi-refractory)

BACT3

lease

nmin
, pmin
B
B

mgC−1

p_qpcPBA

mmolP

p_qlnc, p_qlpc

mmolN mgC−1 ,
mmolP mgC−1

hnB , hBp

p_chn, p_chp

mmolN mgC−1 ,
mmolP mgC−1

βB

p_pu_ea_R3

-

p_version

(-)

Switch

for

bacteria

parameterization:

1=BACT1; 2=BACT2; 3=BACT3

Table 2.4.: Mathematical and code symbols, units and description of the bacterioplankton parameters
(namelist Pelagic_Ecology.nml: PelBac_parameters).

32

2.3. Bacterioplankton
Equation Box 2.3 Bacteria equations
dBc
dt
dBn
dt

=
bio

∑

j=1,6

=
bio

∑

j=1,6

−
dB p
dt

=
bio

dBc
dt

upt

( j)

dBc
dt

Rn

( j)

Rc

dBc
−
( j)
dt
Rc

dBc
−
dt
O(3)

upt
( j)

+

Rc

∑

i=3,4

dBn
dt

lys
(1)

−

Rc

∑

k=5,6

upt

dBn
dt

−
N (i)

prd

dBc
dt

rel

−
N (4)

(2.3.1a)

(k)

Zc

Bn dBc
Bc dt

lys
(1)

+

(2.3.1b)

Rc

prd

Bn
dBc
∑
Bc k=5,6 dt

Zc

( j)

upt

R p dBc
∑ ( j) dt
j=1,6 Rc

rsp

(2.3.1c)

(k)

( j)

Rc

+ fBp

upt,rel

dB p
dt

where d0B is a temperature enhanced background
mortality specific rate and dBd is a density dependent mortality specific rate assumed to be dependent on virus infection.

−
N (1)

B p dBc
Bc dt

lys
(1)

Rc

( j)

f

n,p
R( j)

−

Bp
dBc
∑
Bc k=5,6 dt

( j)

( j)

( j)

R p /Rc R p /Rc
,
= min 1,
opt
opt
p
n

prd
(k)

(2.3.1d)

Zc

!

j = 1, 6 f n,p
( j)
R

(2.3.7)
The
uptake
of
the
nutrient
components
in the
2.3.4. BACT1 parameterization
dissolved and particulate fractions is then derived
(1)
In this parameterization labile DOC (Rc,n,p ) is the from the actual nutrient ratios in the organic matonly class of dissolved organic matter resolved by ter as:
the model and, together with Particulate Organic
(6)
Detritus (Rc,n,p ), constitute the organic substrate
( j)
R dBc upt
dBi upt
available to bacteria.
, i = n, p; j = 1, 6
= i( j)
dt R( j) Rc dt R(c j)
i
(2.3.8)
2.3.4.1. Substrate uptake
The total carbon uptake rate of organic substrate 2.3.4.2. Nutrient release and uptake
in (2.3.1a) is regulated by environmental factors
and substrate availability with a threshold formu- The mineralizers/phytoplankton competitors behavior is controlled by the non-dimensional faclation:
tors fBp and fBn and by the specific relaxation rates
p
n
dBc upt
(2.3.6) νB and νB towards the optimal internal quota. The
∑ dt ( j) = min fBn,p fBT r0B Bc ,
process it involves only phosphate and ammoRc
j=1,6
! nium components:


(1)
n,p
R( j) n,p ( j)
R(1)
Rc
R
+
1
−
f
f
ν
ν
∑ B ( j) c
0B
(1)
j=1,6

R

R

where r0B is the maximum potential specific
( j)
growth rate, νBR is the specific quality-dependent
(1)
uptake rate for substrate, ν0RB is the specific,
quality-independent uptake rate for semi-labile
organic matter and f n,p
is the non-dimensional
R( j)
regulating factor based on the of organic substrate

dB p
dt

dBn
dt

upt,rel

N (1)
upt,rel
N (4)

= fBp νBp

Bp
opt
− pB Bc (2.3.9)
Bc

= fBn νBn

Bn
opt
− nB Bc (2.3.10)
Bc

dBn
dt

upt

=0

(2.3.11)

N (3)

33

2. The pelagic plankton model
In the case of phosphorus as in eq. (2.3.9), for Remineralization dynamics occur only when ininstance, if there is an excess of nutrients in the ternal quotas are in excess with respect to the “opopt
B
cell, Bcp − pB > 0, the non-dimensional parame- timal” value and the dissolved nutrient uptake is
opt
B
ter fBp = −1, and if Bcp − pB < 0 there is direct linearly dependent on the membrane affinity as
uptake from the water as a function of the nu- follows:
trient concentration in a Michaelis-Menten form, dB upt
p
N (1)
=
(2.3.14)
. The same applies to nitrogen,
fBp = N (1)
+hBp
dt N (1)
"
!#
and the direct uptake is regulated by the factor
upt
dB
N (4)
n
p
fB = N (4) +hn .
min a1B N (1) Bc , max 0, popt
B GB −∑
( j)
B
j=1,6 dt R p
2.3.4.3. Excretion
In this parameterization the release term appearing in eq. (2.3.1a) is nil.

dB p
dt

rel

= max 0,
N (1)

∑

j=1,6

dB p
dt

upt
( j)
Rp

− popt
GB
P

!

(2.3.15)

2.3.5. BACT2 parameterization
In this parameterization originally proposed by
Vichi et al. (2004), the substrate available for up(1,6)
(2)
take is composed of Rc,n,p and Rc .

dBn
dt

rel

= max 0,
N (4)

∑

j=1,6

dBn
dt

upt
( j)
Rn

− nopt
B GB

!

(2.3.16)

2.3.5.1. Substrate uptake



dBn upt
hnB
(3)
(4)
a4B n
N + a4B N
Bn ,
(1,6)
∑ dt =min
uptake is parameterized as in BACT1 i=3,4
The Rc
hB + N (4)
N (i)
(2)
!!
(see eq.2.3.6) and the Rc uptake is added.
dBn upt
opt
max 0, nB GB − ∑
( j)
upt
j=1,6 dt Rn
dBc
n,p T
= min fB fB r0B Bc ,
(2.3.12)
(2.3.17)
∑
( j)
j=1,2,6 dt Rc
!


2.3.6. BACT3 parameterization
(1)
n,p
R( j) n,p ( j)
R(1)
R(2) (2)
∑νB fR( j) Rc + ν0B 1 − fR(1) Rc + νB Rc
j=1,6
2.3.6.1. Substrate uptake
(2)

(2)

where νBR is the Rc specific uptake rate. The
uptake of the nutrient components in the dissolved
and particulate fractions is the same defined for
BACT1 in eq. (2.3.8).
2.3.5.2. Nutrient release and uptake

The organic substrate uptake differs from the
BACT1 and BACT2 parameterizations not only
for the wider resolution of the DOM lability/refractivity characteristics given in the introduction,
but also because the DOC uptake is not constrained by the nutritional content of bacterial
cells (see eq. 2.3.3)

The uptake dynamics is similar to the one for phytoplankton, with affinity constants for phosphate
!
and ammonia (a1B and a4B respectively) and also
(
j)
dBc upt
(
j)
∑ dt ( j) = min fBT fBor0B Bc, ∑ νBR Rc
nitrate with a inhibition factor due to ammonium.
Rc
j=1,2,3,6
j=1,2,3,6
This flux regulates the direct uptake and the activ(2.3.18)
ity uptake based on the optimal quota and the net
The uptake of the nutrient components in the
carbon uptake
dissolved and particulate fractions is instead the
dBc upt
dBc rsp
same defined for BACT1 (eq. 2.3.8).
GB = ∑
−
.
(2.3.13)
( j)
dt
dt
(3)
Rc
O
j=1,2,6
34

2.4. Zooplankton
2.3.6.2. Nutrient release and uptake
The Nutrient release and uptake parameterization
is the same described for the BACT1 version (see
Sec. 2.3.4.2)
2.3.6.3. Excretion
This version also defines a bacterial first-order excretion of carbohydrates


B p /Bc
Bn /Bc
= max 0, 1 − opt , 1 − opt
νB Bc
(2)
pB
nB
Rc
(2.3.19)
based on the optimal nutrient content and the constant relaxation time scale νB . The release of semi
refractory DOC is assumed to be at constant rate
(d −1 ), assumed to be proportional to the activity
respiration rates:
dBc
dt

rel

dBc
dt

rel
(3)

Rc

= βB

∑

j=1,2,3,6

dBc
dt

upt
( j)

(2.3.20)

Rc

where βB is the constant fraction of renewal of
capsular material, usually equivalent to about 1/4
of the respiration rate (Stoderegger and Herndl,
1998). Note that this parameter must be set by
the user and a warning is issued if the resulting
bacterial growth efficiency is too low.

2.4. Zooplankton

(6)

• heterotrophic nanoflagellates Zi , protozoa
with dimensions between 2 and 20 µ m,
mainly grazing upon picophytoplankton and
bacteria.
Mesozooplankton is operationally defined in the
model as any zooplankter between 200 µ m and
3 to 4 cm long as an adult, also embracing many species that are traditionally considered part of the microzooplankton when in juveniles stages (Broekhuizen et al., 1995). The
BFM core contains a specific parameterization for microzooplankton and mesozooplankton, but the formal description is presented
for the generic zooplankton as in Vichi et al.
(2007b). The zooplankton parameterization is
modified from Baretta-Bekker et al. (1995) and
Broekhuizen et al. (1995) and it includes the processes of growth due to ingestion and the loss
terms due to excretion/egestion, mortality, respiration and predation. Each zooplankton group
comprises 3 constituents for C, N and P content
as shown in Eq. Box 2.4, however the model
is constructed to reduce the parameterizations to
C-based dynamics assuming constant nutrient to
carbon ratios (see Sec. 7.4).
Si and chl are currently not included as constituents for zooplankton, because biogenic silica
in the form of frustules is directly egested by zooplankters and chl is a negligible part of C and N
in the total biomass of preys.

2.4.1. Food availability

Zooplankton is subdivided into microzooplankton The total amount of food available to zooplankand mesozooplankton with 2 sub-groups defined ton is computed considering the set of possio
n
for each component:
( j)
( j)
as the vector Fi =
ble preys Xi ∈ Pi , Bi , Zi
(3)
• carnivorous mesozooplankton Zi ;
∑X δZ,X eZ,X Xi , where δZ,X is the availability of prey
Xi for predator Z and eZ,X is the capture efficiency.
(4)
• omnivorous mesozooplankton Zi , mainly The product of the latter terms gives the total prefcomprising calanoid copepods ;
erence. There are many definitions of preferences
in the literature, and we have used concepts from
(5)
• microzooplankton Zi , representing the Gentleman et al. (2003) and Gibson et al. (2005)
biomass concentration of microzooplankton to combine the parameterizations described in
with a ESD in the range 20-200 µ m, exclud- Baretta-Bekker et al. (1995) for microzooplanking flagellates and naupliar/larval stages of ton and in Broekhuizen et al. (1995) for mesozoomulticellular zooplankton or meroplanktonic plankton. Availability represents the suitability
larvae of benthic organisms;

35

2. The pelagic plankton model

Equation Box 2.4 Zooplankton equations.
dZc
dt
dZn
dt
dZ p
dt

=
bio

∑

X=P,Z

=
bio

=
bio

dZc
dt

prd

−
Xc

dZc
Fn
∑
Fc X=P,Z dt
Fp
dZc
∑
Fc X=P,Z dt

Symbol

Code

Units

Q10Z

p_q10

(-)

∑

j=1,6

dZc
dt

prd

−
Xc

Xc

( j)

−

Rc

rel

j=1,6

∑

dZn
dt

rel

∑

dZ p
dt

prd

−

rel

j=1,6

dZc
dt

( j)

−

Rn

( j)

Rp

−

rsp

∑

−
O(3)

k=4,5,6

dZn
dt

rel

dZ p
dt

rel

−
N (4)

−
N (1)

dZc
dt

prd

(2.4.1a)

(k)

Zc

Zn
dZc
∑
Zc k=4,5,6 dt

prd
(k)

Zp
dZc
∑
Zc k=4,5,6 dt

prd
(k)

Description
Characteristic Q10 coefficient

p_chuc

mg C

µZ
r0Z
δZ,B
δZ,P
δZ,Z
bZ
ηZ
βZ
εZc
εZn
εZp
nopt
, popt
Z
Z

p_minfood

mg C m−3

p_sum

d−1

Potential specific growth rate

p_paPBA(z,p)

(-)

Availability of pelagic bacteria B to zooplankton Z

p_paPPY(z,p)

(-)

Availability of phytoplankton P to zooplankton Z

p_paMIZ(z,p)
p_srs

(-)
d−1

Availability of microzooplankton P to zooplankton Z
Basal specific respiration rate

p_pu

(-)

Assimilation efficency

p_pu_ea

(-)

Excreted fraction of uptake

p_pe_R1c

(-)

Partition between dissolved and particulate excretion of C

p_pe_R1n

(-)

Partition between dissolved and particulate excretion of N

p_pe_R1p

(-)

Partition between dissolved and particulate excretion of P

p_qncMIZ,

mmolN mgC−1 ,

νZ
d0Z
dZo

(2.4.1c)

Zc

hF
Z

(2.4.1b)

Zc

m−3

Michaelis constant for total food ingestion
Feeding threshold

Maximum nutrient quota

mgC−1

p_qpcMIZ

mmolP

1

d−1

Specific rate of nutrients and carbon excretion

p_sd

d−1

Specific mortality rate

p_sdo

d−1

Oxygen-dependent specific mortality rate

p_chro

mmolO2 m−3

Half saturation value for oxygen

Table 2.5.: Mathematical and code symbols, units and description of the microzooplankton parameters
(namelist Pelagic_Ecology.nml: MicroZoo_parameters).

36

2.4. Zooplankton
Symbol

Code

Units

Description

Q10Z

p_q10

(-)

Characteristic Q10 coefficient
Potential specific growth rate
Specific search volume

r0Z

p_sum

d−1

vZ

p_vum

m3 mg C −1 d−1

δZ,P
δZ,Z
δZ,Z
bZ
ηZ
βZ
nopt
, popt
Z
Z

p_paPPY(z,p)

(-)

Availability of phytoplankton P to zooplankton Z

p_paMIZ(z,p)

(-)

Availability of microzooplankton P to zooplankton Z

p_paMEZ(z,p)
p_srs

(-)
d−1

Availability of mesozooplankton P to zooplankton Z
Basal specific respiration rate

νZ
d0Z
dZdns
γZ

p_puI

(-)

Assimilation efficency

p_peI

(-)

Excreted fraction of uptake (faeces production)

p_qncMEZ,

mmolN mgC−1 ,

Maximum nutrient quota

p_qpcMEZ

mmolP mgC−1

1

d−1

Specific rate of nutrients and carbon excretion

p_sd

d−1

Specific mortality rate

p_sdo

m3 mgC−1 d−1

Density-dependent specific mortality rate

p_sds

(-)

Exponent for density dependent mortality

p_clO2o

mmolO2 m−3

Half saturation value for oxygen

Table 2.6.: Mathematical and code symbols, units and description of the mesozooplankton parameters
(namelist Pelagic_Ecology.nml: MesoZoo_parameters).
r

of the prey and is assumed to be mostly dependent on the prey’s nominal dimensions. Capture
efficiency (or relative preference) is also a nondimensional factor which is set to 1 for mesozooplankton and it is density-dependent in microzooplankton, eZ,X = XcX+cµZ , according to the threshold
half-saturation density µZ (µZ = 0 for mesozooplankton).

plankton (hFZ = v0Z ), because this parameter is
Z
generally available in the literature. For brevity,
in the zooplankton equations we will use the following notation to indicate the total ingestion rate
in units of each constituent:

2.4.2. Ingestion

2.4.3. Excretion/egestion

The first term on the right hand side of eq.
(2.4.1a) is the total carbon ingestion, which corresponds to the sum of all the predation loss
terms in the carbon equations of the other functional groups preyed by zooplankton. Applying
the inter-functional group conversion defined in
eq. (1.0.2), the rate term for each predation processes is parameterized with a Type 2 formulation
(Gentleman et al., 2003),

Metabolic rates in zooplankton are assumed to
be closely coupled to growth, therefore total ingested carbon is used part for net production,
part for respiration and the remainder is egested/excreted. The parameters that can be measured in laboratory experiments are net growth
efficiency ηZ and the egested portion of ingested
material βZ . The ingestion rate in eq. (2.4.3)
is not directly affected by prey quality in our
present formulation, although this process appears to modulate zooplankton feeding rates substantially (cf. Mitra and Flynn, 2005). Nevertheless, the definition of constant (optimal) nutrient quota in zooplankton (Baretta-Bekker et al.,

prd

prd

δZ,X eZ,X Xc Fc
Zc
Fc
Fc + hFZ
Xc
Zc
(2.4.2)
which is traditionally rewritten in terms of the
specific search volume in the case of mesozoodZc
dt

=−

dXc
dt

= fZT rZ0

Ij = ∑
X

dZ j
dt

prd

j = c, n, p.

(2.4.3)

Xj

37

2. The pelagic plankton model
1997), equivalent to the Threshold Elemental Ratios of Andersen et al. (2004, TER), implies that
the ingestion of low-quality (i.e. nutrient-poor)
food lead to the disposal of the ingested carbon in
excess, thus effectively limiting biomass growth.
On the other hand, an excess of nutrients, as
for instance due to the ingestion of phytoplankton under “luxury uptake” conditions, leads to an
increase of the nutrient remineralization rates as
shown below in eqs. (2.4.14-2.4.15). The release
of extra C is parameterized as an increase of the
egestion rates of organic carbon compounds or,
alternatively, by increasing the respiration rates.
Both processes are well documented in freshwater
zooplankton (Frost et al., 2004; Anderson, 2005)
and we have decided to parameterize increased
excretion rates. The two pathways are equivalent
from the point of view of internal element regulation in zooplankton, but the consequences of one
choice or another on the biogeochemical cycling
of carbon are still to be investigated both experimentally and in model studies.
The carbon loss term in (2.4.1a) thus represents
the sum of the activity excretion/egestion (higher
for mesozooplankton because of sloppy feeding),
the mortality rates and the nutrient-limited excretion of organic carbon. For the microzooplankton
the carbon loos term reads:

∑

j=1,6

dZc
dt

rel
( j)
Rc

= βZ Ic + (d0Z + dZo (1 − fZo )) fZT Zc

(2.4.4)
for mesozooplankton the carbon loss term
reads:

density valid only for mesozooplankton (dZdns = 0
for microzooplankton).
The balancing flows of C, N, P, Qc,n,p
are comZ
puted from the actual elemental ratios of ingested
material:
ΓiZ =

(1 − βZ ) Ii
,
ηZ Ic

i = n, p

(2.4.6)

the (1 − βZ ) factor is the assimilation fraction of
element uptake. The SWITCH (see Tab. 2.7)
indicates which is the limiting element, default
is carbon, if the ΓiZ is lower than the internal
quota Zi /Zc indicates that i is limiting. The ratios between the ΓiZ and the respective zooplankton internal elemental quota of Zi /Zc are crosscompared, the lowest defines the most limiting element i. For example in the case of nitrogen limitation (N column Tab.2) the correction on carbon
(QcZ ) is the difference between effective ingestion
of carbon (ηZ Ic ) and the effective ingestion of ni):
trogen scaled by the nitrogen TER (nopt
Z
QcZ = ηZ Ic −

(1 − βZ )
In ,
nopt
Z

(2.4.7)

2.4.4. Respiration
Taking into account the energy cost of ingestion
(1 - assimilation - egestion) and basal metabolism
the total respiration rate can be written as:
dZc
dt

rsp
O(3)

= (1 − ηZ − βZ )Ic + bZ fZT Zc

(2.4.8)

where the constant basal respiration rate bZ is also
considered. The energy cost of ingestion is comdZc rel
γ
o
T
dnsputed
c
considering
the following metabolic bal= βZ Ic +(d0Z +doxy0Z (1− fZ )) fZ Zc +dZ ZcZ +Q
Z
dt R(6)
ance:
c
(2.4.5)
The released fraction is further divided into
Ic = Gc + Ec + Rc
(2.4.9)
particulate (faecal pellets) and dissolved organic
forms using a constant percentage εZc (mesozoo- where Ingestion (I) is partitioned in Growth (G),
plankton is assumed to have no dissolved prod- Excretion (E) and Respiration (R). Assimilation
ucts). Mortality is parameterized as senescence efficiency ηZ is:
with a first-order rate based on a constant d0Z ,
Gc
(2.4.10)
ηZ =
as oxygen regulated component doxy0Z , and as
Ic
a grazing closure by higher trophic levels not resolved in the model, which is a power function of and excretion is:

38

2.4. Zooplankton
Limiting
Element ->

C

SWITCH

default

QcZ

0

QnZ

(1 − βZ ) In − nopt
ηZ Ic
Z

0

ηZ (Ic − QcZ )
(1 − βZ ) In − nopt
Z

QZp

(1 − βZ ) I p − popt
ηZ Ic
Z

(1 − βZ ) Ip − popt
ηZ (Ic − QcZ )
Z

0

N

if

ΓnZ
ZN /ZC

<

ΓZp
Z p /ZC

ηZ Ic −

P

and

(1−βZ )
nopt
Z

ΓnZ
ZN /ZC

<1

if

ΓZp
Z p /ZC

<

ΓnZ
Zn /ZC

ηZ Ic −

In

and

(1−βZ )
popt
Z

ΓZp
Z p /ZC

<1

Ip

Table 2.7.: Mesozooplankton formulation to eliminate the excess of the non-limiting constituent. The
SWITCH control determines whether C, P or N is the limiting element. ηZ is the net
and popt
are
growth efficiency factor, βZ is the egested portion of ingested material, nopt
Z
Z
the optimal Threshold Elemental Ratios (TERs).

2.4.6. Inorganic nutrients
Ec = βZ Ic

(2.4.11)

The third terms on the right hand side of eqs.
(2.4.1b) and (2.4.1c) parameterize the zooplankton excretion of inorganic nutrients, which occur
(2.4.12) only when the internal nutrient quota exceed the
Rc = (1 − ηZ − βZ )Ic
opt
optimal quota for P and N, popt
Z and nZ , respectively. The following formulation is applied to
2.4.5. Excretion/egestion of organic
microzooplankton:
nutrients
after some algebra we can express R in term of I:

The nutrient dynamics for zooplankton given in
eqs. (2.4.1b) and (2.4.1c) are mainly derived from
carbon dynamics taking into account the nutrient content of the total food uptake. The excretion/egestion rate of organic nutrients is obtained
from eq. (2.4.4):



Zp
opt
max 0,
− pZ
Z(2.4.14)
p
Zc
N (1)


dZn rel
Zn
n
opt
=ν max 0,
− nZ
Z(2.4.15)
n
dt N (4) Z
Zc

dZ p
dt

rel

=νZp

and the time scales of excretion are controlled by
the specific constant rates νZp and νZn . The excre
tion is in the form of phosphate and urea, but the
dZi rel
Zi 
T
dns γZ
β
=
I
+
d
f
Z
+
d
Z
∑ dt ( j) Zc Z c 0Z Z c Z c
latter in the model is assumed to be as labile as
Ri
j=1,6
the ammonium, therefore the rate is directed to
i = n, p
(2.4.13) the N (4) pool. In the case of mesozooplankton the
for microzooplankton is subsequently partitioned inorganic excretion is congruent with the formubetween particulate and dissolved matter accord- lation for carbon excretion (Tab.2):
ing to the non-dimensional fraction εZi , which parameterizes the different distribution of nutrients
dZ p rel
= d0Z fZo fZT Z p + QZp
(2.4.16)
between structural parts and cytoplasm. Mesodt N (1)
zooplankton only releases particulate organic detritus.

39

2. The pelagic plankton model
losses due to pelagic chemical reactions:
dZn
dt

rel
N (4)

= d0Z fZo fZT Zn + QnZ

(2.4.17)

2.5. Non-living components

dO(2)
dt

( j) gpp

3

= Ωoc ∑

j=1

bio

−Ωoc fBo

2.5.1. Oxygen and anoxic processes
The dynamics of dissolved oxygen and carbon
dioxide are important closures of global biogeochemical cycles. We do not describe here the
exchange of gases at the air sea interface which
is assumed to be a purely physical process and
has been thoroughly investigated elsewhere, especially for CO2 (Olsen et al., 2005).
Anaerobic processes and denitrification dynamics are a consequence of oxygen dynamics
and are described here for completeness, although
they are of limited impact in the well-oxygenated
euphotic zones of the open ocean. Nevertheless,
these processes are important for the sulfur cycle
and for the fate of exported carbon in the mesoand bathypelagic layers of the ocean, where bacteria are the major drivers of these processes. To
account for hypoxic and anoxic remineralization
in the water, the original ERSEM parameterization of anaerobic processes in the sediments proposed by Ruardij and Van Raaphorst (1995) was
extended to the pelagic system by Vichi et al.
(2004). The state variable “reduction equivalents”
N (6) (Table 1.1 and Fig. 1.1) is an inorganic
state variable containing all the reduced chemical species and assumed to be chemically equivalent to the sulphide ion HS− . The basic constituent is indicated with the letter R because this
variable account for all the reduced biochemical
products, although it should be mostly regarded
as sulphur S. Reduction equivalents are produced
as a result of bacterial anoxic respiration and are
partly used for the parameterization of denitrification processes and partly for direct sulphide production.
The pelagic net production of oxygen is derived
from the sum of gross primary production and
community respiration rates from phytoplankton,
zooplankton and bacteria, also subtracting the

40

dPc
dt
dBc
dt

6

−Ωoc ∑

j=4

O
rsp

O(3)

!

dt

+

O(3)

nit

reox

N

sinkr

o

1 dN (6)
− r
Ωo dt
(3)

(2.5.1)

All the rates are converted into oxygen units
by means of constant stoichiometric coefficients.
Since bacteria are active both under aerobic and
anaerobic conditions the bacterial oxygen demand (eq. 2.3.4) is partitioned into oxygen consumption and reduction equivalent production by
using the oxygen regulating factor fBo in (2.3.2).
The nitrification rate is a source term of the nitrate
equation (2.5.6b), and a sink term for ammonium
(eq. 2.5.6c) and oxygen (eq.2.5.1). Nitrification
is not explicitly resolved but parameterized with a
simple first-order dependence on ammonium and
oxygen concentrations:
dN (4)
dt

nit

= Λnit(4) fnT
N (3)

N

O(2)
N (4)
O(2) + ho

(2.5.2)

where Λnit
is the constant specific nitrification rate
N4
and fnT a temperature regulating factor with the
Q10 formulation shown in (eq. 2.1.1).
The formation of reduction equivalents is parameterized converting the biological oxygen demand of bacteria (under low oxygen conditions)
into sulphide ions by using the stoichiometric coefficient Ωro as:
dN (6)
dt

=Ωro Ωoc 1 − fB1o
bio


 dBc
dt

dN (3)
− Ωo Ωn
dt
r

eo

+

+

O(3)
( j) rsp
dZc

dN (4)
−Ωn
dt

( j) rsp

dPc
−
dt
(3)

denit

sinkn

rsp

+
O(3)

dN (6)
−
dt

reox

sinkr

(2.5.3)

2.5. Non-living components
The utilization of nitrate as an electron acceptor in microbial metabolic reactions is parameterized in an indirect way. Firstly, when the oxygen
level falls below the threshold level and fBo < 1
(eq. 2.3.2), the metabolic formation of reduction
equivalents begins according to the carbon mineralization rate (eq. 2.3.4). The denitrification reaction is favored with respect to the strictly anaerobic sulphate reduction, therefore a portion of this
oxygen demand is redirected towards the denitrification process. In order to achieve this net effect, the changes in the redox conditions enhance
the denitrification flux in the following way:
dN (3)
dt

denit
denit

=Λ

N (3)

T

fn




 dBc
1 o
Ωc 1 − fBo
∗
Mo
dt

rsp
O(3)



N (3) ·

(2.5.4)

where Λdenit
is the specific denitrification rate at a
N (3)
reference anoxic mineralization Mo∗ . If nitrate is
still present in the water, the bacterial rate of production of reduction equivalents N (6) is converted
to nitrate consumption, mimicking the bacteriamediated denitrification reactions. This chemical
rate leads to a direct production of gaseous N2 in
the water, which is the only time rate of change
for state variable O(4) in the model.
Furthermore, as long as there is some oxygen
left, reduction equivalents are also quickly reoxidized at the following rate:
dN (6)
dt

reox

= Λreox
(6)
N

sinkr

O(2)
N (6)
(2)
O + ho

(2.5.5)

where Λreox
is the (constant) specific daily reN (6)
oxidation rate, and ho is the half-saturation oxygen concentration. When oxygen and nitrate are
completely depleted the last two terms in (2.5.3)
become zero and the process turns to a strict
anaerobic formation of sulphide ions coupled to
the availability of the organic substrate.

ganic CFFs for dissolved compounds are considered here (Fig. 1.1): phosphate, nitrate (nitrate
+ nitrite), ammonium and silicate with the equations shown in Box 2.5 and derived from the processes described in the previous sections.
The pelagic cycle of phosphate N (1) in (eq.
2.5.6a) is affected by phytoplankton uptake (eq.
2.2.23), bacterial uptake/release (eq. 2.3.9) and
excretion from zooplankton groups (eq. 2.4.14).
The pelagic processes for nitrate N (3) shown
in eq. (2.5.6b), involve phytoplankton uptake
described in eq. (2.2.1b) and the nitrification
and denitrification process parameterizations described in equations (2.5.2) and (2.5.4), respectively.
Ammonium (eq. 2.5.6c) is consumed by phytoplankton as described in eq. (2.2.23) and remineralized (or utilized) by bacteria according to
the quality of the substrate and their internal content of nitrogen according to eq. (2.3.10). Zooplankton participates in the ammonium dynamics
through the excretion of urea, which is assumed
to be directly available in the form of ammonium,
as shown in eq. (2.4.15).
The pelagic cycle of silicate is quite simple
in the model because of the many uncertainties
linked to the complex dynamics of this element
in the water. Silicate concentration was originally
only affected by diatom uptake (eq. 2.2.1e), but
a simple first-order reaction parameterizing bacterial dissolution (e.g. Bidle and Azam, 2001) has
been introduced accounting for the dissolution of
silicate frustules as:
(6) rmn

dRs
dt

(6)

= Λrmn
f T(6) Rs
s
N (5)

R

(2.5.7)

where Λrmn
is the constant specific dissolution
s
T
rate and f (6) is the temperature regulating factor
R
as in eq. (2.1.1), mimicking bacterial activity enhancement at higher temperatures.
Iron is made available in dissolved form
through remineralization of biogenic particles
2.5.2. Dissolved inorganic nutrients
produced by phytoplankton and zooplankton. The
The pelagic cycles of dissolved inorganic nutri- biochemical pathways of the remineralization
ents are essential components of any biogeochem- process are not completely clear and involve
ical model of the marine ecosystem. Five inor- both siderophores and photochemical reactions.

41

2. The pelagic plankton model

Symbol
Ωoc
Ωon
e on
Ω
Ωro
Ωrn
Λnit

Code

Units

Description
mgC−1

MW_C

mmolO2

p_qon_nitri

mmolO2 mmolN−1

p_qon_dentri

mmolO2 mmolN−1
−

Unit conversion factor and stoichiometric coefficient
Stoichiometric coefficient for nitrification reaction
Stoichiometric coefficient for denitrification reaction

p_qro

mmolHS

p_qro * p_qon_dentri

mmolHS − mmolN−1

Stoichiometric coefficient nitrogen-anoxic reaction

p_sN4N3

d−1

Specific nitrification rate at 10o C

Q10n
ho

p_q10N4N3
p_clO2o

(-)
mmolO2 m−3

Q10 factor for nitrification/denitrification reaction
Half saturation for chemical processes

hr

p_clN6r

mmolHS − m−3

Half saturation oxygen concentration for anoxic pro-

Λdenit

p_sN3O4n

d−1

Specific denitrification rate

Mo∗

p_rPAo

mmol O2 m−3 d−1

Reference anoxic mineralization rate

Λreox
N6

p_rOS

d−1

Specific reoxidation rate of reduction equivalents

Q10N5

p_q10R6N5

(-)

Q10 factor for dissolution of biogenic silica

Λrmn
s

p_sR6N5

d−1

Specific dissolution rate of biogenic silica

εPAR
λW
cR(6)

p_PAR

(-)

Fraction of Photosynthetically Available Radiation

p_eps0

m−1

Background attenuation coefficient

p_epsR6

m2 mg C−1

C-specific attenuation coefficient of particulate detritus

vR6

p_sediR6

m d−1

Settling velocity of particulate detritus

N4

mmolO2

−1

Stoichiometric coefficient for oxic-anoxic reaction

cesses
N3

sed

Table 2.8.: Chemical stoichiometric coefficients and general parameters involving pelagic
components.

Equation Box 2.5 Dissolved inorganic nutrient equations.
dN (1)
dt
dN (3)
dt
dN (4)
dt
dN (5)
dt

42

4

= −∑
bio

j=1
4

= −∑
bio

j=1
4

= −∑
bio

( j) upt

dPp
dt

p

+ fB

N (1)
( j) upt
dPn

dN (3)
+
dt
(3)

dt

N
( j) upt
dPn

dt

j=1

bio

N

p

+ fB

N (4)

(1) upt

dPs
= −
dt

dB p
dt

(k) rel

upt,rel

+
N (1)

nit

N

dBn
dt

∑

k=4,5,6

dN (3)
−
dt
(4)

(2.5.6a)
N (1)

denit

(2.5.6b)
sinkn
(k) rel

upt,rel

+
N (4)

dZ p
dt

∑

k=4,5,6

dZn
dt

N

dN (4)
−
dt
(4)

nit

(2.5.6c)
N3

(6) rmn

dRs
+
dt
(5)

(2.5.6d)
N (5)

2.5. Non-living components
Equation Box 2.6 Dissolved organic matter equations
(1)

dRc
dt

3

=

dt

j=1

bio

(1)
dRi

3

=

=
bio

∑

j=1

bio

(2)

dRc
dt

∑

dPc
dt

( j) exu

dPc
dt

dBc
−
dt
(1)

Rc
( j) exu
dPi

dt
exu

(1)

−

Ri

(1)

Rc

(1)

Ri

dBc
−
(2)
dt
Rc

upt
(2)

Rc

(1)

dBc
dt

dZc
dt

∑

+

Rc

k=5,6

(1)

(2.5.8a)
(1)
Rc

(k)

upt

+

Rc

dBc
+
dt

Since all these processes are primarily bacterialmediated, it is assumed here that dissolved Fe is
released from detritus according to a first-order
relationship as for silicate (2.5.7):
(6)

(k) rel

upt

(k)
k=5,6 Zc

rel
(2)

Rc

∑

Zi

!

(k) rel

dZc
dt

i = n, p

(2.5.8b)

(1)

Rc

(2.5.8c)

1
−1 and
with a given time constant Λscv
f = 40 years
with the further assumption that scavenging results in adsorption onto sinking particles and consequently sequestration in the deeper layers.

rmn

dR f

(6)

= Λrmn
f T(6) R f
f

dt

R

(2.5.9) 2.5.3. Dissolved and particulate

N (5)

where Λrmn
is a constant specific dissolution rate
f
T
and f (6) is the temperature dependence. Both
R
numbers are currently unknown, and therefore
they need to be adjusted by trial-and-error for balancing the iron cycle in the ocean. The inclusion
of iron as an explicit component of zooplankton
and bacteria may link this process to the direct
excretion of organisms and bacterial regeneration
activity, once the important pathways and timescales have been properly assessed by laboratory
and in situ experiments.
Dissolved inorganic iron species are scavenged
onto particle surfaces owing to hydroxide precipitation. Since the concentration of iron ligands
is about 0.6 nM in the deep ocean, Johnson et al.
(1997) suggested that iron scavenging can be parameterized with a constant rate when the iron
is above this threshold. Ligand dynamics have
been further investigated by Archer and Johnson
(2000); Parekh et al. (2004); Lefevre and Watson
(1999), but the simplest approach as proposed
by Johnson et al. (1997) and Aumont et al. (2003)
has been used here:
scv


dN (7)
(7)
= Λscv
min
0,
N
−
0.6
(2.5.10)
f
dt

organic matter
The state variables describing dissolved organic
matter shown in Box 2.6 have been introduced in
Sec. 2.3 because they are tightly linked to the parameterizations of pelagic bacteria. Three biogeochemical basic constituents C, N and P and are
thus described by 3 equations shown in. DOM is
produced by phytoplankton, bacteria and microzooplankton and used as organic substrate by bacteria. The different degrees of lability of DOM are
(1)
reflected in the nutrient content of R j , which regulates bacterial uptake as shown in eq. (2.3.18),
(2)
(3)
while variables Rc and Rc only contains the
carbon constituent. The equation (2.5.8c) for
semi-labile DOC (carbohydrates) is derived from
the production terms implemented in the different
bacteria parameterizations (Secs. 2.3.4, 2.3.5 and
2.3.6). Bacteria are allowed to release carbohydrates only in BACT2 and BACT3 parameterizations.
Particulate detritus is instead described by 4
equations, one for each biogeochemical basic
constituent C, N, P, Si as in Box 2.7. The carbon,
nitrogen and phosphorus component of particulate detritus in equations (2.5.11a) and (2.5.11b),

sink f

43

2. The pelagic plankton model
respectively, are produced by all the members of −log10 ([H + ])), dissolved inorganic carbon conthe planktonic community except bacteria, which centration (DIC) and total alkalinity (TA). These
are the only utilizers of this component according species are governed by the following relations:
to eq. (2.3.18).
+
[HCO−
3 ] · [H ]
The pelagic cycle of biogenic silica is in(2.5.13)
K1 =
[CO2 ]
stead restricted to the release of diatom frustules
through mortality and other lysis processes as in
+
[CO2−
3 ] · [H ]
eq. (2.2.28) and via micro/mesozooplankton pre(2.5.14)
K2 =
[HCO−
3]
dation (including sloppy feeding) with the addition of the chemical dissolution shown in eq.
(2.5.7).
2−
DIC = [CO2 ] + [HCO−
(2.5.15)
3 ] + [CO3 ]
Particulate iron dynamics are the consequence
of processes described in equations (2.2.9),
[CO2 ]
(2.5.9) and (2.5.10). Particulate organic Fe is also
pCO2 =
(2.5.16)
derived from zooplankton egestion and mortalK0
ity. It is assumed that zooplankton is never ironlimited and the iron fraction of the ingested phyto2−
−
TA = [HCO−
3 ] + 2[CO3 ] + [B(OH)4 ] +
plankton is directly egested as particulate detritus.
3−
+[OH −] + [HPO2−
]+
4 ] + 2[PO4 (2.5.17)

2.5.4. The carbonate system

+
−
+[H3 SiO−
4 ] − [H ]F − [HSO4 ] +

The aquatic chemistry of inorganic carbon forms
(state variable O(3) ) is a further extension to the
original ERSEM formulation and it is activated
with the key INCLUDE_PELCO2. The theory
of dissolved inorganic carbon chemical reactions
is well understood (Zeebe and Wolf-Gladrow,
2001) and we propose here a slightly revised formulation with respect to Blackford and Burkill
(2002), which also takes into account the algorithms proposed by the Ocean Carbon Model Intercomparison Project (Doney et al., 2004).
In the ocean, inorganic carbon is present
in three different forms: free carbon dioxide
([CO2 ] = [CO2 ]aq + [H2CO3 ]), bicarbonate ion
2−
(HCO−
3 ) and carbonate ion (CO3 ). These carbonate species reach the following equilibrium:

−[HF] − [H3 PO4 ]

+
2−
+
CO2 + H2 O ⇋ HCO−
3 + H ⇋ CO3 + 2H
(2.5.12)
defined by the equilibrium constants K1 and K2
for the first and second reaction, respectively. The
carbonate system in seawater is described in terms
of 7 chemical species, i.e., free carbon dioxide
(CO2 ), bicarbonate ion (HCO−
3 ), carbonate ion
(CO2−
),
carbon
dioxide
partial
pressure
in seawa3
ter (pCO2 ), hydrogen ion concentration (pH =

44

(2.5.18)

The species appearing in this latter equation are
expressed in terms of their equilibrium constants
and their total elemental concentrations. Total alkalinity is therefore computed as a function of:
TA = f ([H + ], DIC, K1 , K2 , Kw, Kb ,
K1p , K2p , K3p , Ksi , Ks , K f , bt, st, f t,
(2.5.19)
O(3) , O(5) , N (5) , N (1) )

(2.5.20)

where K1 and K2 are the previously seen equilibrium constants for carbonic acid (H2CO3 ) and bicarbonate ion (HCO−
3 ), K0 is the Henry’s constant which regulates CO2 solubility in seawater,
Kw is the ion product of water, Kb is the dissociation constant for boric acid (B(OH)3 ), K1p , K2p
and K3p are the dissociation constants for phosphoric acid (H3 PO4 ), di-hydrogen phosphate ion
2−
(H2 PO−
4 ) and hydrogen phosphate ion (HPO4 )
respectively, Ksi is the dissociation constant for
silicic acid (Si(OH)4 ), Ks is the dissociation constant for bisulfate ion (HSO−
4 ), K f is the dissociation constant for hydrogen fluoride (HF),

2.5. Non-living components
Equation Box 2.7 Particulate organic detritus equations.
(6)

dRc
dt

4

=
bio

(6)
dRi

dt

4

=

∑

dBc
−
dt
(6)

Rc
( j) lys
dPi

dt

j=1

bio

=
bio

dPs
dt

−

Ri

(6)

Rc

(6)

(1) 6

+
(6)
Rs

(6)

Ri

(1) lys

(6)

dRs
dt

∑

j=1

( j) lys

dPc
dt

Ps

(1)

Pc

∑

j=4

upt
(6)

Rc

(k) rel

6

+∑

dBc
dt

k=4

upt

6

(6)

Rc

+∑

(1) prd

dPc
dt

dZc
dt

(2.5.11a)
(6)
Rc

(k)

Zi

(k)
k=4 Zc

i = n, p

(2.5.11b)

(6)

Rc

(6) rmn

dRs
−
dt
( j)

Zc

(k) rel

dZc
dt

(2.5.11c)
N (5)

bt is the total boron concentration ([B(OH)3 +
[B(OH)−
4 ]), st is the total sulfate concentration
−
([HSO4 ] + [SO2−
4 ]), f t is the total fluoride concentration ([HF] + [F − ]), and O(3) , O(5) , N (1)
and N (5) are the model state variables for DIC,
total alkalinity, total phosphorus concentration
2−
3−
([H3 PO4 ] + [H2 PO−
4 ] + [HPO4 ] + [PO4 ]), and
total silica concentration ([Si(OH)4 ]+ [H3 SiO−
4 ]),
respectively. Note that total fluoride, sulfate and
boron concentrations ( f t, st and bt) are given as
dependent on the temperature and salinity only.
Currently, the model does not take into account
changes of alkalinity due to calcification.
The effects of biogeochemical processes on alkalinity linked to the N cycle are taken into account with a correction to the total alkalinity value
after Wolf-Gladrow et al. (2007). The uptake of
1 mole of N leads to (i) an increase of alkalinity
by 1 mole when nitrate is the source, (ii) to a decrease of alkalinity by 1 mole when ammonia is
used. Nitrification leads to a decrease of TA by
2 moles per mole of nitrate formed and denitrification leads to an increase of TA by 1 mole per
mole of nitrate converted. Thus, considering the
equations in (2.5.6b) and (2.5.6c), the source form
equation for TA reads:
dO(5)
dt

4

=
bio

( j)

∑

dPn
dt

−2

dN (3)
dt

j=1

The equilibrium constants are also computed
from temperature and salinity according to the
methods described in Zeebe and Wolf-Gladrow
(2001) where the total pH scale is used. The
carbonate equilibrium calculation is done in units
of mol · kg−1 , therefore conversions from model
units are done by taking into account the actual
seawater density.
This is a five equation system (eqs. 2.5.132.5.17) with seven unknown variables: the system is therefore determined when two of the seven
variables are known, which in this case are DIC
and TA. Total alkalinity equation is used to compute [H + ] by means of the Newton-Raphson convergence method. A dedicated diagnostic variable
is also introduced to store the total hydrogen ions
and to allow the restart of the model from previous simulations. The other carbon species are
eventually computed from the carbonate equilibrium equations.
Alternatively, the model allows the usage of a
simplified computation of DIC equilibria as suggested by Follows et al. (2006), which does not
imply an iterative procedure. This solution is suggested for coupled configurations because in that
case the inaccuracy of the computation is less important due to the presence of advective and difupt
( j) upt
4
fusive processes.
dPn
−∑
+
Finally, the biological production and conj=1 dt
N (3)
N (4)
sumption of CO2 considered in the model can be
nit
denit
easily derived by collecting the first 4 terms on the
dN (3)
+
(2.5.21) right hand side of eq. (2.5.1) without considering
dt
(4)
N
sinkn
the stoichiometric factor Ωoc and taking the total

45

2. The pelagic plankton model
bacterial respiration as

controlled with a specific namelist found in file
Carbonate_Dynamics.nml and listed in Tab. 2.9.
3
dO(3)
The namelist also allows to read an external file
−
+
=∑
dt
dt
dt
j=1
(3)
(3)
containing data for atmospheric concentration of
bio
O
O
(k) rsp
rsp
inorganic carbon or to set it constant to a certain
dZc
dBc
− ∑
−
value of the mixing ratio. Atmospheric partial
dt O(3) k=4,5,6 dt
O(3)
pressure can also be computed or read from an
(2.5.22) external file.
The parameters of the carbonate system are
( j) gpp
dPc

46

( j) rsp
dPc

!

2.5. Non-living components

Parameter

Units

Description

pCO2Method

Integer Switch

Methods for computing pCO2
1 = AtmCO20* sea level pressure (slp0 in Param.nml
or external input AtmSLP_N),
2 = Compute the pCO2 using the Magnus formula:
Sea Level Pressure and Dew Point Temperature have
to be provided

AtmCO20

ppmv

Set constant atmospheric CO2 concentration (ppm),
also referred as the CO2 mixing ratio

calcAtmpCO2

Logical

Compute the atmospheric pCO2

phstart

pH

Initial guess of the seawater pH.

K1K2

Integer Switch

Acidity constants parameterization
1- Roy et al. (1993); DOE (1994); pH on total scale;
2- Mehrbach et al. (1973) refit by Dickson & Millero
(1987); pH on Sea Water Scale;
3 - Mehrbach et al.(1973) refit by Lueker et al.
(2000), pH on total scale;
4 - Hansson (1973b) data as refitted by Dickson and
Millero (1987), pH on Sea Water Scale.

MethodCalcCO2

Integer Switch

Method for Carbonate System computation: 1 Approximate static solution; 2 - Standard OCMIP
iteration; 3 - Approximation by Follows et al. (2006).

CalcBioAlkFlag

Logical

Compute the corrections of alkalinity due to

M2XACC

Threshold

Accuracy of the iterative scheme for OCMIP (with

M2PHDELT

pH

pH range for the root search in the OCMIP scheme

M2MAXIT

Integer

maximum number of iterations for OCMIP (with

biogeochemical processes
MethodCalcCO2=2)
(with MethodCalcCO2=2)
MethodCalcCO2=2)
Caconc0

mol

m−3

Calcium ion concentration (In "Seawater : Its
composition, properties and behaviour", Open
University Course Team, 1995)

Canorm

Logical

Normalize Calcium ion concentration by sea water
salinity

AtmCO2_N

ppmv

Read external file with atmospheric CO2 mixing
ratio data.

AtmSLP_N

Pa

AtmTDP_N

oC

Read external file with Sea Level Pressure data (with
pCO2Method=1 and 2).
Read external file with atmospheric Dew Point
Temperature data (with pCO2Method = 2).

Table 2.9.: Parameters controlling the carbonate system (from namelist Carbonate_Dynamics.nml)

47

2. The pelagic plankton model

48

3. The sea ice biogeochemical model
3.1. A model for sea ice
biogeochemistry
Sea ice is a rich habitat for microbial community.
The most abundant species found are unicellular
microalgae, mostly diatoms. When sea ice forms,
many organisms are either passively or actively
entrapped in the salty brines. Their rate of survival in the new habitat depends on their adaptation and/or acclimation to the new environmental conditions (low temperature, high salinity and
low light intensities) and on the external supply of
nutrients and gases from seawater. Some organisms may die in isolated brines, some may survive or encyst, some may find a favorable habitat and actively grow. Concentrations up to 1000
mg m−3 of diatom chlorophyll have been found
in Antarctic sea ice (Thomas and Dieckmann,
2002).
The sea ice biogeochemical cycle is also
strongly related to its oceanic counterpart. This is
of extreme biological importance at the end of the
ice season, when sea ice starts melting and a sea
ice algae bloom occurs. The fate of this biomass
depends on the rate of melting and on the vertical
stability of the water column. If the stratification
is high and the rate of melting is low, sea ice algae
may stay long time in the upper part of the water column and may seed a pelagic phytoplankton
bloom. Polar blooms represent a relevant fraction
of the carbon production in some regions of the
world, such as the Ross Sea and the Weddell Sea
in Antarctica, and the Barents Sea in the Arctic.
If the rate of melting is high and the stratification
is low, the sea ice biomass may rapidly sink to the
bottom of the ocean and likely become a sink for
the atmospheric CO2 . In both cases, the size and
weight of the organisms affect the sinking velocity.
Sea ice biota has been studied for only a few
decades. Few regions have been highly character-

ized, but sea ice biological variability at different
temporal and spatial scales is still lacking. Sea
ice is one of the largest ecosystem on earth, but
is also one of the less sampled: sampling sea ice
biota is in fact not an easy task. It is costly and
time-consuming and often it is done in severe environmental conditions. In absence of data and
remote sensing facilities, modelling can help understanding the sea ice ecosystem, as well as it
can provide the wider picture of its qualitative and
quantitative importance, which is still missing.
The sea ice extension of the BFM described here (BFM-SI) has been published in
a series of papers (Tedesco et al., 2010, 2012;
Tedesco and Vichi, 2014) and this chapter is a
modified version of Tedesco and Vichi (2010).
The reader interested in a full theoretical description and a review of the existing techniques for
sea ice biogeochemical modelling is addressed to
Tedesco and Vichi (2014). This latter paper also
describes the default example available as part of
the BFM presets, STANDALONE_SEAICE.

3.2. Model structure
This implementation focuses on primary producers, which are the most abundant group of organisms found in sea ice and the most relevant group
in terms of export of biomass to the ocean. The
biogeochemical equations of the sea ice algae dynamics are written according to the formulation
explained in Chapter 1.
The model layout (Fig. 3.1) takes advantage of
the same biological processes of the pelagic BFM.
The focus is here on primary producers, which
are assumed to differently adapt to the new physical environment. The main differences between
BFM and BFM-SI stand in the type and number
of functional groups (Table 3.1), in the parameters
assigned to several physiological and ecological

49

3. The sea ice biogeochemical model
processes (Table 3.2) and in the dimensional form tative of sea ice primary producers:
they represent. While pelagic state variables are
• Adapted diatoms, which are meant to be
expressed in terms of their constituent per cubic
highly adapted to the environment and also
meters, the BFM-SI state variables are expressed
show distinct skills in acclimation. They are
in terms of constituent per square meters. The
supposed to be first light-limited and, later
strategy of coupling will be further described.
in the bloom, dependent on nutrient availBFM-SI totally resolves 28 state variables (Fig.
ability. They have an Equivalent Spherical
3.1, Table 3.1):
Diameter (ESD) of 20-200 µ m and preyed
• 2 sea ice algae functional groups (adapted
by adult mesozooplankton (> 200 µ m) and
diatoms and surviving sea ice algae, mostly
microzooplankton of larger dimensions (20represented by autotrophic nanoflagellates)
200 µ m), which are not currently present
in
the sea ice system, but act in the pelagic
• 6 inorganic variables for nutrients and gases
BFM when sea ice melts and algae are re(phosphate, nitrate, ammonium, silicate,
leased in the water column. Sea ice diatoms
oxygen and carbon dioxide)
are the main source of biogenic silica and
• 2 organic non-living groups for dissolved
differ from the other subgroup being their
and particulate detritus
growth limited by dissolved silicate.
• 1 generic group of aerobic and anaerobic sea
ice bacteria
• 1 generic group of sea ice fauna
Each state variable interacts with the others
through the universal physiological and ecological processes depicted in Fig. 3.1. This component already includes parameterizations of sea
ice bacteria and fauna, which follow the same dynamics as their counterpart in the pelagic realm.
Nonetheless, both groups are technically available for future studies.
As for the pelagic model, nitrate is assumed
here to be the sum of both nitrite and nitrate.
All the nutrient:carbon ratios in chemical organic
and living functional groups are allowed to vary
within their given range and each component has
a distinct biological time rate of change. This
kind of parametrization is meant to mimic the
adaptation of organisms to the diverse availability of nutrients and light observed in nature, and
also allow to recycle organic matter depending on
the actual nutrient content.

• Surviving sea ice algae, which may be
mostly represented by autotrophic nanoflagellates, are meant to only survive in the sea
ice environment, being less adapted to it and
showing lower skills of acclimation. However, they may be able to grow in sea ice if
the diatoms bloom is quickly exhausted - for
instance, for depletion of silicate - and a sufficient amount of nutrients is still available
for their growth. Their ESD is 2-20 µ m and
are mainly externally preyed by pelagic microzooplankton.

The mathematical notation used here is the same
defined for the pelagic BFM (Chap 1). Sea ice
algae are involved in several processes: gross primary production (gpp), respiration (rsp), exudation (exu), cell lysis (lys), nutrient uptake (upt),
predation (prd) and biochemical synthesis (syn).
Both subgroups share the same form of primitive equations, but are differentiated in terms of
the values of the physiological parameters (Table
3.2). There are 5 living CFFs that describe the
constituents of the generic variable sea ice algae
A (with constituents C, N, P, Si and Chl, see Ta3.3. Sea ice Algae Dynamics
ble 3.1) and thus for each group we have 4 or 5
As a first implementation of the BFM in sea ice, 2 equations as shown in Equation Box 3.1:
distinct subgroups have been chosen as represen-

50

3.3. Sea ice Algae Dynamics

Figure 3.1.: Scheme of the state variables of BFM-SI and interactions within BFM-SI and with external systems.
The rate of change of carbon in sea ice algae depends on gross primary production, exudation, respiration, lysis and predation (Eq. 3.3.1a).
Gross primary production is the rate of change
of sea ice algae carbon Ac due to photosynthesis, which involves an uptake of dissolved carbon
dioxide F (3) . It is written as:
dAc
dt

fAT = QA10

gpp

= fAT fAE fAs rA0 Ac

tiplicative, non-dimensional regulating factors for
temperature, light and silicate, which vary from 0
to 1.
Temperature is regulating several physiological processes. Its effect is expressed in a nondimensional form by fAT :
T −10
10

(3.3.3)

(3.3.2)

where QA10 is the characteristic doubling temperature parameter (Table 3.2).
where the rA0 is the maximum specific photosynMany relevant biological processes, such as pothetic rate under nutrient-replete, light saturated
tential
photosynthesis, are also affected by the
conditions (Table 3.2). The f functions are mulF (3)

51

3. The sea ice biogeochemical model

Table 3.1.: List of the sea ice model state variables. Legend: IO = Inorganic; LO = Living organic;
NO = Non-living organic. The subscript i indicates the basic components of the group, e.g.
(1)
(1)
(1)
(1)
(1)
(1)
Ai ≡ (Ac ; An ; A p ; As ; Al ).
Variable

Type

Components

Description

I (1)

IO

P

Phosphate (mmol P m−2 )

I (3)

IO

N

Nitrate (mmol N m−2 )

I (4)

IO

N

Ammonium (mmol N m−2 )

I (5)

IO

Si

Silicate (mmol Si m−2 )

F (2)

IO

O

Dissolved oxygen (mg C m−2 )

F (3)

IO

C

Carbon dioxide (mg C m−2 )

(1)

LO

C N P Si Chl

Adapted diatoms (mg C m−2 , mmol N-P-Si m−2 , mg Chl-a m−2 )

Ai

(2)

LO

C N P Chl

Surviving sea ice algae (mg C m−2 , mmol N-P m−2 ,mg Chl-a m−2 )

Ti

LO

CNP

Sea ice bacteria (mg C m−2 , mmol N-P m−2 )

Xi

LO

CNP

Sea ice fauna (mg C m−2 , mmol N-P m−2 )

(1)

NO

CNP

Dissolved organic detritus (mg C m−2 , mmol N-P m−2 )

(6)

NO

C N P Si

Particulate organic detritus (mg C m−2 , mmol N-P-Si m−2 )

Ai

Ui

Ui

non-dimensional light regulating factor fAE :
λbio , where:
!
2
(6)
( j)
EPAR
E
λbio = ∑ cA Al + cU (6) Uc .
fA = 1 − exp −
(3.3.4)
EK
j
where EPAR is the Photosynthetic Available Radiation (PAR). EPAR is parametrized according
to the Lambert-Beer formulation with depthdependent extinction coefficients:
EPAR(z) = εPAR Fsw e(λs +λi )z+

R0
z

λbio (z′ )dz′

(3.3.5)

Thus, λbio takes into consideration the extinction
due to sea ice algae chlorophyll and to particulate detritus, while dissolved substances and other
inorganic matter are not currently taken into account. The cA and cU constants are the specific
absorption coefficients of each suspended substance (Table 3.2).
EK is the light saturation parameter, that is the
ratio between the maximum chl-a specific photosynthetic rate and the maximum light utilization
coefficient, i.e.:

where Fsw is the short-wave surface irradiance
flux and may be derived from data or from
a coupled physical model, such as the one of
Tedesco et al. (2009, 2010). The irradiance flux
is then converted by BFM-SI from W m−2 to the
P∗
units of µ E m−2 s−1 as done in the pelagic model.
EK = m∗ .
α
εPAR is the coefficient determining the portion of
PAR in Fsw . Light propagation takes into account As for pelagic phytoplankton of BFM:
the extinction due to the background extinction of
Ac
snow/sea ice λs,i and due to particles in the sea ice
Pm∗ = fAT fAs rA0
Al

52

(3.3.6)

(3.3.7)

(3.3.8)

3.3. Sea ice Algae Dynamics

Table 3.2.: Ecological and physiological parameters in BFM-SI.
Symbol

A(1)

A(2)

Description

r0A

1.5

2.0

Maximum specific photosynthetic rate (d−1 )

QA10

2.0

2.0

Characteristic Q10 coefficient (-)

0
θchl

0.035

0.03

Optimal quotum chl-a:C (mg chl mg C−1 )

0
αchl

1.8 e−3

3.8 e−6

Maximum light utilization coefficient (mg C (mg chl)−1 mE−1 m2 s)

ds

0.1

-

Half saturation value for Si-limitation (mmol Si m−2 )

bA

0.05

0.1

Basal specific respiration rate (d−1 )

γA

0.10

0.10

Activity respiration fraction (-)

βA

0.05

0.20

Excreted fraction of primary production (-)

dAp,n,s

0.1

0.2

Nutrient stress threshold (-)

dOA

0.1

0.1

Maximum specific lysis rate (d−1 )

a1

2.5 10−3

2.5 10−3

Specific affinity constant for P (m−2 mg C−1 d−1 )

a3

2.5 10−3

2.5 10−3

Specific affinity constant for N-NO3 (m−2 mg C−1 d−1 )

a4

2.5 10−3

2.5 10−3

Specific affinity constant for N-NH4 (m−2 mg C−1 d−1 )

sopt
A(1)

0.03

-

Standard Si:C ratio in sea ice diatoms (mmol Si mg C−1 )

smax
A(1)

0.085

-

Maximum Si:C ratio in sea ice diatoms (mmol Si mg C−1 )

pmin
A

1.97 10−4

1.97 10−4

Minimum phosphorus quota (mmol P mgC−1 )

popt
A

7.86 10−4

7.86 10−4

Optimal phosphorus quota (mmol P mgC−1 )

pmax
A

1.57 10−3

1.5710−3

Maximum phosphorus quota (mmol P mgC−1 )

nmin
A

3.78 10−4

3.78 10−4

Minimum nitrogen quota (mmol N mgC−1 )

nopt
A

1.26 10−3

1.26 10−3

Optimal nitrogen quota (mmol N mgC−1 )

nmax
A

2.52 10−3

2.52 10−3

Maximum nitrogen quota (mmol N mgC−1 )

cA

10.0 e−3

10.0−3

Specific absorption coefficient for chlorophyll-a (m2 (mg Chl-a)−1 )

cU

0.1 e−3

0.1e−3

Specific absorption coefficient for detritus (m2 (mg C)−1 )

Ωoc

1
12

1
12

Unit converison factor and stochiometric coefficient (mmol O2 mgC)−1

53

3. The sea ice biogeochemical model
Equation Box 3.1 Sea ice algae equations
dAc
dAc
= = − ∑
dt
j=1,6 dt
dAn
dt
dA p
dt

=

=

dAl
dt

dAc
( j) dt

Uc
upt

∑

dA p
dt

−
I ( j)

=

F

dAc
−
dt
(3)

dAn
j=1,6 dt

∑

lys

I

Up

dA p
− ∑
j=1,6 dt
(1)

dAs
dt

= θchl

gpp

upt

(1) upt

(1)

dAs
dt

dAn
j=3,4 dt

lys

I

resp

Uc

F (3)

dAc
−
dt
(1)

(3.3.1a)

lys

(3.3.1b)
( j)

Un

(3.3.1c)
( j)

(1) lys

dAs
−
dt
(5)

dAc
dt

exu

(3.3.1d)
(6)
Us

gpp

exu

F

(1)
Uc

dAc
−
dt
(3)

!

dAc
−
dt

rsp

lys

F

(6)
Uc

dAc
+
dt
(3)

!

Al
Ac

(3.3.1e)

(3.3.9) where n(p)opt
A is the nitrate(phosphate) optimal
min
where fAT is the regulating factor for temper- ratio, while n(p)A is the nitrate(phosphate) minature, fAs is the regulating factor for silicate, imum quota (Table 3.2).
The respiration rate of Eq. (3.3.1a) is written
rA0 is the maximum specific photosynthetic rate
as:
under nutrient-replete, light-saturated conditions
0
α ∗ = fAT fAs αchl

0 is the maximum slope of the productionrsp
and αchl
dAc
= fAT bA Ac +
irradiance curve at optimal conditions (Table 3.2).
dt (3)
s
The fA is parametrized as an external limiting
F
gpp
exu !
factor with a Michaelis-Menten form:
dAc
dAc
(3.3.13)
−
γA
I (5)
dt (3)
dt (1)
s
(3.3.10)
fA = (5)
F
Uc
I + ds
where ds is the Michaelis-Menten constant for where fAT is the metabolic regulating factor for
temperature, bA is a constant specific rate of resSiO2 uptake inhibition (Table 3.2).
piration and γA is a fraction of the assimilated proThe exudation rate of Eq. (3.3.1a) reads:
exu
duction (Table 3.2).
dAc
The loss of carbon via lysis of Eq. (3.3.1a) is
= [βA +
dt (1)
written
as:
Uc

(1 − βA )(1 −

dAc
fAn,p )]
dt

gpp

(3.3.11)
F (3)

dAc
∑ dt
j=1,6

lys

=
( j)

Uc

1
dO Ac (3.3.14)
fAp,n + dAp,n A

where βA is a constant fraction of carbon uptake
(Table 3.2) and fAn,p is a Liebig-like regulating where dAp,n is the nutrient stress threshold and dOA
factor for internal nutrient ratio:
is the maximum specific lysis rate (Table 3.2).
min
A
/A
−
n
n
c
A
,
fAn,p = min
min
nopt
−
n
A
A
!
A p /Ac − pmin
A
(3.3.12)
min
popt
A − pA
54

3.4. Nutrient Supply and Dynamics
The chlorophyll rate of change of Eq. (3.3.1e)
is due to chlorophyll synthesis. The net chlorophyll synthesis is a function of acclimation to
light conditions, nutrient availability and turnover
rate. As in BFM, it is assumed that nutrientstressed cells releasing substantial amounts of dissolved organic carbon tend to regulate their internal chl:C ratio in order to avoid unconstrained decreases.
The rate of change of net photosynthesis is thus
primarily controlled by the dynamical chl:C ratio
θchl proposed by Geider et al. (1998), which regulates the amount of chl-a in the cell according to
a non-dimensional ratio between the realized photosynthetic rate in Eq. (3.3.1e) and the maximum
potential photosynthesis, i.e.:

in BFM-SI as a smaller value for the half saturation of silica and a larger value for the standard
Si:C quotum in adapted diatoms (Table 3.2).
The boundary fluxes are currently added as additional source terms to the biogeochemical equations and solved explicitly. For instance, in the
case of an inorganic nutrient in sea ice (e.g. nitrate, I (3) ), the complete equation is written as
(1) upt

dI (3) dAn
=
dt
dt

I

(2) upt

dAn
+
dt
(3)

I

dI (3)
+
dt
(3)

f lux

(3.4.1)
N (3)

where the first two terms on the right hand side
represent the uptake from sea ice algae and the
last one is the flux of nutrient at the boundaries.
The external mechanisms of nutrients replenishment (exchange with the ocean and atmospheric
fAE rA0 Ac
0
(3.3.15) deposition) will be analyzed in the next sections.
θchl = θchl 0
αchl EPARAl
The uptake of nitrogen and phosphorous by algae
(Eq. 3.4.1) is regulated by a Droop kinetics:
0 is the maximum quotum chl-a :C and
where θchl
0 is the maximum slope of the productionupt
αchl

dAn
a3A I (3)
=
min
irradiance curve at optimal growth conditions (Ta∑ dt
i=3,4
I (i)
ble 3.2). The same considerations about down
regulation and chlorophyll losses as detailed in
+a4A I (4) Ac , nopt
A GA
Sec. 2.2.6 for phytoplankton are valid for sea ice
!

An 
algae.
max
T 0
+ fA rA nA −
(3.4.2)
Ac
Ac

3.4. Nutrient Supply and
Dynamics

Nutrients supply for algal growth comes from the
mixed layer up to the ice sheet for sustaining bottom communities, but also from snow deposition
through brine drainage for surface communities
and from in situ regeneration processes.
Even in isolated brine pockets, bacteria,
heterotrophic protozoa and small metazoans
have been shown to regenerate the major
nutrients (Arrigo et al., 1995), but not silicate.
Silicate dissolution and regeneration
may be slower than demand and can be
the major limiting factor for diatoms growth
(Lizotte and Sullivan, 1991), shifting the community from being diatom-dominated to flagellatesdominated (Dieckmann et al., 1991). The slow
regeneration of silicate in sea ice is parametrized

dA p
dt

upt

= min a1A I (1) Ac , popt
A GA
I (1)

+ fAT rA0



Ap 
Ac
pmax
A −
Ac

!

(3.4.3)

where the a constants are the membrane affinity
for nitrate, ammonium and phosphate (Table 3.2).
The uptake of silicate is, instead, only function
and of the net proof the maximum Si:C ratio smax
A(1)
duction GA(1) of Eq. (3.3.1a):
dAs
dt

upt

= smax
G
A(1) A(1)

(3.4.4)

I (5)

Whenever sea ice algae carbon is lost by lysis,
a proportional loss is found for algae nutrient content and is distributed between a dissolved and a

55

3. The sea ice biogeochemical model
(6)
particulate fraction. For instance, the equations of
Particulate organic matter (POM, Ui in BFMthe lysis rate for phosphorous are:
SI) is made of C, N, P and Si:
lys

dA p
dt

dA p
dt

lys
(1)
Up

=

∂ Ac
∑ ∂t
j=1,6

(6)

∂ Uc
∂t

(3.4.5)

∂t

(6)

Up

Ap
=
Ac

lys

∂ Ac
pmin
A

(6)

Uc

lys

lys

Uc

(6)
Up

∂ Ap
−
∂t
( j)

(6)

.(3.4.6)

∂ Ui
∂t

bio

∂ Ac
=∑
j=1 ∂ t

2

bio

( j) lys

2

(3.5.3)
(6)
Uc

( j) lys

∂ Ai
=∑
j=1 ∂ t

i = n, p, s (3.5.4)
(6)

Ui

Silicate is instead only released in particulate where the silicate component of POM is only
valid for the release of sea ice diatoms frustules.
form:
(1) lys

∂ As
dt

(6)

Us

(1)

(1) lys

As ∂ Ac
= (1)
∂t
Ac

.

(3.4.7)

3.5. Gases and Detritus
The following equations are derived by combining terms from the previous sections in order to
ensure mass conservation. The net production of
oxygen is due to the gross primary production and
to algae respiration rates:
!
( j) gpp
j rsp
2
∂ F (2)
∂
A
∂
A
c
c
= Ωoc ∑
−
∂t
∂
t
∂
t
j=1
(3)
(3)
bio

F

F

Ωoc

is the stoichiometric conversion factor
where
to oxygen units in respiration and photosynthesis
(Table 3.2).
(1)
Dissolved organic matter (DOM, Ui in BFMSI) is a non-living functional group including C,
N and P constituents. In this current implementation, DOM is produced by sea ice algae (Eq.
3.3.11), though in the complete setup shown in
Fig. 3.1 , DOM will be produced also by sea ice
bacteria and microzooplankton:
(1)
∂ Uc

∂t

(1)

∂ Ui
∂t

56

2

=
bio

2

bio

∑

j=1

( j) exu
∂ Ac

∂t

(3.5.1)

(1)
Uc

( j) exu

∂ Ai
=∑
j=1 ∂ t

i = n, p.
(1)

Ui

3.6. The coupling strategy

(6)

Uc

(3.5.2)

A full description of the coupling issues is given
in Tedesco and Vichi (2014); here we briefly
mention some major aspects. The construction of
a biological system in sea ice implies the coupling
with the underlying biology of the ocean. BFMSI is in fact coupled to a simplified version of the
pelagic lower trophic levels in a surface layer of
ice-covered oceans. The standard pelagic BFM
has been simplified in a way that every sea ice
group has its own pelagic counterpart and there is
no loss of material between the two systems (Fig.
3.2). The number of functional groups is reduced
and the total number of state variables computed
are 34. The included groups and subgroups in the
pelagic BFM are:
• 2 phytoplankton (diatoms and autotrophic
nanoflagellates)
• 3 zooplankton (omnivorous mesozooplankton, microzooplankton and heterotrophic
nanoflagellates)
• 1 pelagic bacteria
• 9 inorganic variables for nutrients and
gases (phosphate, nitrate, ammonium, silicate,oxygen, carbon dioxide)
• 2 organic non-living variables for dissolved
and particulate detritus.

3.6. The coupling strategy

Figure 3.2.: Structure of the current coupling between physics and biogeochemistry (boxes with continuos line). The sea ice physical model passes offline the Environmental Sea Ice Variables (ESIV) to BFM-SI. BFM-SI is online coupled to BFM through the exchange of
Non Living Inorganic (NLI), Non Living Organic (NLO) and Living Organic (LO) matter. Boxes with dashed lines represent the possibility for the sea ice physical model to
be also coupled to an ocean physical model, which would pass the Environmental Ocean
Variables (EOV) to BFM.
Differently to BFM-SI, the pelagic zooplankton includes 3 different subgroups and mesozooplankton may effectively control the fate of the
sea ice algae released into the water and the
magnitude of the phytoplankton bloom. Consequently, the diversity of feeding behaviors of zooplankton is maintained.
The exchange of matter between the ocean system and the sea ice system will be the subject
of the next section. We focus here on the coupling strategy. BFM-SI has been built as a layer
model, in a similar way as the benthic component of BFM is designed. A layer model means
that the model is two-dimensional (concentrations
are expressed in units of mass per square meters):
all the members of a generic layer of the system,
as for instance the benthic oxic layer, have to be
thought as their constituents are homogeneously
distributed on a stratum at the interface with the
ocean. However, the state variables of the pelagic
BFM are expressed as concentrations in a volume of sea water. The thickness of the considered layer is taken into account in order to have
every CFF of both systems (pelagic ad benthic) in

the same units per cubic meters. This concept can
be applied also to the coupling between BFM and
BFM-SI, considering that the sea ice biology is
distributed in a certain layer of sea ice (the sea ice
bio, BAL of Fig. 3.3). The thickness of the layer
is computed by a physical sea ice model, such as
the one developed by Tedesco et al. (2010), which
is a time-varying layer whose thickness depends
on the physical properties of sea ice (temperature,
salinity and brine volume).

3.6.1. Boundary fluxes: the sea
ice-ocean interface
The major exchanges between the sea ice and the
pelagic systems are the fluxes of organic and inorganic matter at the interface. Ice structure determines the porosity and therefore the rate of exchange: frazil ice is less porous than congelation
ice and algae may experience, for instance, nutrient depletion in the former.
The rates of advection and diffusion of dissolved and particulate substances from seawater
into the porous bottom layer of ice may also depend on under-ice current velocities, tides and

57

3. The sea ice biogeochemical model

Figure 3.3.: Model representation of the sea ice and pelagic coupling during a generic ice season. The
pelagic BFM exchange matter and gases with the SEA ICE BIO layer.
atmospheric pressure cycles (Cota et al., 1991).
Additionally, high algal biomass in the platelet
layer has been found to reduce the flux of nutrients to algae communities of the upper congelation layer in Antarctica Lizotte and Sullivan
(1991).
High biomass is usually found in area of slow
ice growth rate (Legendre et al., 1992). During
the growth season, convection in the skeletal layer
enhances nutrient fluxes (Cota et al., 1991). During the melting season, the supply of freshwater
increases the stratification just below the ice and
reduces the flux of nutrients upward by reducing mixing and friction velocity (Gosselin et al.,
1985). Even if the enlargement and interconnections of the brine channels allows a greater
biomass accumulation and nutrient supply, it also
leads to an increase in the biological loss by cells
sinking in the water column.
It is assumed here that the entrapment of dissolved matter follows the same partitioning of
salt in sea ice, that is the dynamics of dissolved
constituents is treated as the salinity dynamics,
as proposed by Tedesco et al. (2009), considering
the concentration of dissolved matter in seawater and sea ice growth rate. Fluxes are defined
positive upward, i.e. towards sea ice in Fig. 3.2.
The nitrate flux is described here as an example of
an inorganic nutrient exchange. The flux is composed of a positive part (entrapment) during sea
ice formation when sea ice growth is positive, and

of a negative part during the melting phase:
dI (3)
dt

f lux−oce

N (3)

∂ hi
= max 0,
∂t



·

!
I (3)
max 0, N −
hbio

 (3)
∂ hi I
+ min 0,
∂ t hbio
(3)

(3.6.1)

where N (3) is the nitrate concentration in seawater, hi is the ice thickness and hbio is the thickness
of the biologically active layer in sea ice. When
ice starts melting (the last term in the equation
above), the release to the water column depends
only on the sea ice melting rate and on the sea ice
concentration. The same flux, converted into units
of volume concentration, is included with opposite sign in the dynamical equation for pelagic nutrients.
The entrapment of particulate matter is assumed to be only a function of the seawater concentration, the sea ice growth rate and the actual
available space in the sea ice matrix (brine volume, Vbio ), while the release during the melting
phase is parametrized as function of the sea ice
melting rate and sea ice concentrations, as in Eq.
(3.6.1). For instance, the flux of the chlorophyll
(1)
component of sea ice algae Al , from the pelagic
(1)
variable Pl to sea ice is defined as:
(1) f lux−oce

dAl
dt

58



P(1)




∂ hi
(1)
= max 0,
Pl Vbio
∂t
 (1)

∂ hi Al
. (3.6.2)
+ min 0,
∂ t hbio

3.6. The coupling strategy

3.6.2. Boundary fluxes: the sea
ice-atmosphere interface

munities, while during snow flushing episodes
the accumulated nutrients in the melting snow
may become available also for bottom communities, if brines are permeable. This latter process is the one that it is currently considered in
the model for the parametrization of the flux of
atmospheric nutrients to bottom sea ice, as shown
in Fig. 3.2. Phosphate and nitrogen concentration in snow may be derived from an atmospheric
model or may be prescribed using available observations. The flux of nutrients due to precipitation
becomes then a linear function of the snow melting rate, similarly to the flux of matter from sea
ice to the ocean during melting of sea ice, i.e. for
nitrate:

Nutrient content in snow and liquid precipitation
may directly lead to additional nutrient availability for surface communities, and indirectly to internal and bottom communities, and finally to surface waters once sea ice has totally melted away.
In the Baltic Sea, for instance, 5% of the total
annual flux of N and P and 20–40% of lead and
cadmium are deposited as snow: due to the intense stratification of most of the Baltic waters,
sea ice is the major source of nutrients and trace
elements to the surface Baltic waters during the
melting season (Granskog et al., 2006).
The physical processes responsible of the
f lux−atm


bioavailability of atmospheric nutrients to the sea
∂ hs
dI (3)
=
min
0,
N snow (3.6.3)
ice community are snow ice formation and snow
dt
∂
t
N snow
flushing. When snow ice forms the fraction of
snow that mixes with flooding seawater may bring where N snow is the nitrate concentration in snow
an additional source of nutrients for internal com- and h is the snow thickness.
s

59

3. The sea ice biogeochemical model

60

Part II.

BFM code description

61

4. Installation, configuration and
compilation
The default basic configuration of BFM is the STANDALONE model, which simulates the biogeochemistry of a water volume with given depth. The typical example is a micro- or mesocosm batch
culture. Only the pelagic system can be currently simulated with this configuration. The model is
forced with prescribed or analytical functions of temperature, salinity and light and can be integrated
for any desired period with three different numerical schemes. The model output is in NetCDF.

4.1. Installation
The BFM code repository is made available through a distributed git repository, which requires access
authentication by using the username and password provided in the registration to the BFM website
(http://www.bfm-community.eu/). The code can be downloaded using the following git instruction
% git clone

http://username:password@dev.cmcc.it/git/bfm.git

Alternatively, it is possible to retrieve a tarball archive from the website. After downloading the
tarball file, create a directory for the BFM and move into that directory. Uncompress and extract the
contents of the tarball. The downloaded file will have a different name based on the version. For
example:
%
%
%
%

mkdir $HOME/BFM
cd $HOME/BFM
gunzip bfm-release-.tgz
tar xvf bfm-release-.tar

This operation will create and populate the directory bfm containing the source files and the directory run for running the various test cases which are described in Sec. 5.5. The full directory tree of
the BFM core is given in Fig. 4.1.
The minimal user local settings to use BFM configuration/compilation tools is provided by means
of the environmental shell variable:
• BFMDIR The root directory of the BFM model (the path to the installed bfm- directory).

4.1.1. System Requirements
The BFM can be run on any UNIX-based architecture and it has been tested on linux, Mac OSX and
IBM-iDataplex systems.
The software requirements are:
• PERL (version 5.8.2 and above), used automatically during compilation time to generate the
code from the model template (see Sec. 7).

63

4. Installation, configuration and compilation

Figure 4.1.: Directory tree of the current BFM version.

64

4.2. Configuration: BFM Presets
• FORTRAN 90/95 compiler. Under linux and Mac OSX the model can be currently compiled
with gfortran (version 4.5 and above) and ifort (version 12 and above).
• NetCDF
library
(version
4
and
above
are
recommended,
http://www.unidata.ucar.edu/software/netcdf) . It is mandatory that the library has been
compiled with the same compiler used for the model compilation since the F90 netcdf module is
used.

4.2. Configuration: BFM Presets
An experiment is defined by a model layout structure, namelists, and initial input values. These
definitions put together with compilation and run options is what is called a PRESET. A preset is
contained together in a folder and it is all you need to create a particular experiment.
The basic operations for configuration, compilation, and deployment of the model are performed
automatically by the script bfm_configure.sh in the folder ${BFMDIR}/build/, while the folder containing the presets is ${BFMDIR}/build/configurations/.
Execute the following command to show a list of all available presets:
./bfm_configure.sh -P
Presets are composed of three different main files:

- configuration (compilation and deployment options)
This file uses the F90 namelist format to set values. Below is reported an example of the configuration
file format:
&BFM_conf
MODE
= ’STANDALONE’,
CPPDEFS = ’BFM_STANDALONE INCLUDE_PELCO2 INCLUDE_DIAG’,
ARCH
= ’gfortran.inc’,
PROC
= 1,
EXP
= ’standalone.pelagic’, /
...
=
/
Not all the options available in Tab. 4.1 can be specified in the configuration file and all the values
are overridden by the command line options. Available options are: MODE, CPPDEFS, BFMDIR,
NEMODIR, ARCH, CLEAN, PROC, NETCDF, EXP, NMLDIR, PROC and QUEUE. Do not use the
character " to surround values, but use ’ instead.
There are 2 major macros that belong to the core of the BFM source code and additional compilation
keywords that activates model components. The first 2 macros are “negative”, which means that their
options are activated by default and needs to be deactivated when desired.
Note that, the MODE field is used to specify if the code is compiled for standalone or coupled
execution, thus meaning that specific settings of the script generator will be used. Beside the STANDALONE mode, additional fields will be made available following the release of the official coupling
support for external hydrodynamic models.

65

4. Installation, configuration and compilation
BFM_STANDALONE: This macro controls the main memory structure of the BFM. In default configuration the BFM allocates its own memory array. When the macro is defined at compilation
(-DBFM_STANDALONE) all the variables belonging to the main memory are defined as pointers
and are not allocated during the initialization phase.
BFM_NOPOINTERS: by default, all the F90 symbols of the BFM state variables are defined as pointers to the main memory (cf. Sec. 6.2). By defining this macro, it is possible to substitute all
the pointers with the actual memory location. For instance, all the occurrences of the keyword
P1c(:) which are pointers to the memory address of phytoplankton carbon content, are substituted with the name D3STATE(ppP1c,:). To facilitate this feature, it is requested that all the
references to the pointers in the code are written with the indication of the number of dimensions:
P1c should always be indicated as P1c(:), otherwise a pre-compiler error is generated when
BFM_NOPOINTERS is defined. This macro is strongly recommended when BFM is coupled to
hydrodynamical models.
INCLUDE_PELCO2: Define and activate the pelagic variables of the carbonate system.
INCLUDE_PELFE: Add iron component and activate the related dynamics in the system.
INCLUDE_DIAG: Activate the computation and storage of user-defined diagnostics described in
Chap. 8.
EXPLICIT_SINK: Activate the explicit usage of serapate source and sinks arrays (Sec. 6.2).

- layout (memory layout configuration file)
This file is the template to define the model components layout (Sec. 7) and user-defined diagnostics
(Sec. 8.3). For a detailed description of the BFM memory layout content refer to Sec. 6.2.

- namelists_bfm (namelist file with standard values)
This file contains all the standard values of the namelists1 used for the experiment. Namelists are
checked for consistency against the source code at generation time and the files effectively used for
the simulation will be copied in the $BFMDIR/run/ directory.
During this phase the BFM named constants corresponding to state variables (O2o, N1p, P1c,
...) are substituted by the progressive number which is a function of the chosen layout file. This
check will guarantee an executable file compatible with the namelists provided. The regular user will
generally work with generated namelists and usually there is no need to change any keyword in this
file besides the values in case there is the willingness to generate a new standard set of values.
After running the experiments, it is possible to modify generated namelists in the deployed experiment folder ( $BFMDIR/run/$EXP ). If you want to preserve these changes into a new "preset" you
can use the tool "nmlcreator". This is for instance desired when a new group is added as in Sec.
7.2 and the parameters of the new group have to become default for the preset. To execute this tool do
the following steps
cd $BFMDIR/build/scripts/conf
./nmlcreator.pl -i NAMELIST_IN -o NAMELIST_OUT

where
1 BFM

66

input namelists follows F90 namelist standard

4.3. Compiling the BFM STANDALONE
• NAMELIST_IN: Complete path of the file containing namelist/s to read (e.g.. $BFMDIR/run/standalone.pelagic/Pelagic_Ecology.nml)
• NAMELIST_OUT:
Complete
path
of
the
"namelist"
file
of
the
preset
where
to
append
at
the
end
of
the
input
(e.g.:
$BFMDIR/build/configurations/STANDALONE_PELAGIC_Z/namelists).
Basically, nmlcreator adds the "filename_nml_conf" key with the name of the file at the end
of each namelist read. It is important to note that it will not delete the existing namelists in the target
preset, it will only append the new namelist at the end. In this case, an autimatic back substitution of
variables numerical indexes with the corresponding variable names (e.g. O2o, N1p) is not possible
and it has to be done manually by the user if required.

4.3. Compiling the BFM STANDALONE
Configuration and deployment of the model is done automatically by the script bfm_configure.sh. After setting-up $BFMDIR variable in your shell environment, you can go to $BFMDIR/build directory
and execute bfm_configure.sh using default preset, namely the STANDALONE_PELAGIC:
./bfm_configure.sh -gcd
bfm_configure.sh generates, compiles and/or deploys files needed for the BFM model, as specified
in the preset folder. There are three different steps to have an executable model (see Fig. 4.2):
g

(Generation) generate fortran files (.h , .f90) based on layout information and the single
namelist files (.nml) from namelists_bfm in $BFMDIR/build/tmp/ and
also place fortran files in the main code tree ($BFMDIR/src);

c

(Compilation) compile the fortran code usig files of the main code tree ($BFMDIR/src)
and create the executable file (.x) in $BFMDIR/bin;

d

(Deployment) create experiment folder in $BFMDIR/run/ and copy
namelists and executable for the simulation.

At least one of these three basic instruction is mandatory.
An execution without errors2 will show the output directory for running the experiment and the
command to execute. Output after successful compilation should end with:
.................................. Makefile is ready.
STANDALONE_PELAGIC generation done!
Starting STANDALONE_PELAGIC compilation...
STANDALONE_PELAGIC compilation done!
Go to $BFMDIR/run/standalone.pelagic and execute command: ./bfm_standalone.x
Presets are loaded using option “-p” followed by the preset name, e.g.:
./bfm_configure.sh -gcd -p CO2TEST
For a full description of the script usage and available options, see Tab. 4.1 or execute the following
command:
./bfm_configure.sh -h
2 If

you experience problems during the generation or compilation, visit http://bfm-community.eu/
documentation for help

67

4. Installation, configuration and compilation

Figure 4.2.: Steps for experiment deployment

Option

Argument

VERBOSE
PRESET
MODE
CPPDEFS
BFMDIR
NEMODIR
ARCH
PROC
FAST
NETCDF

-v
-p
-m
-k
-b
-n
-a
-r
-f
-t

EXP
NMLDIR
PROC
QUEUE

-x
-l
-r
-q

Description
Compilation options
Verbose mode to print all messages
Preset to generate the configuration
Mode for compilation and deployment
Key options to configure the model
path to root directory of BFM
path to root directory of NEMO (only for NEMO coupling)
Specify compilation Architecture file
Number of procs used for compilation (and execution)
Fast mode. Don’t execute "clean" command in compilation
Path to netcdf library and header files
Deployment options
Name of the experiment for generation of the output folder
Input dir where are the namelists to run the experiment
Number of procs used for running (same as compilation)
Name of the queue number of procs used for running

Table 4.1.: List of optional arguments for bfm_configure

68

5. Running the STANDALONE model
5.1. Description
The execution of the configuration script described in the previous section creates the executable and the running environment for the default preset of the STANDALONE model (called
STANDALONE_PELAGIC). The list of the available presets is obtained with the command
bfm_configure.sh -P and the results of the test cases are detailed in Sec. 5.5. The generated namelists are copied to the directory $BFMDIR/run/standalone.pelagic and the model
is run by executing ./bfm_standalone.x or ./bfm_standalone.x &> outputfile to
redirect the output messages to a file in bash.
Most of the BFM settings can be controlled via namelists. The STANDALONE model behaviour
is controlled with 2 major namelist files BFM_General.nml and Standalone.nml. The first
file contains the namelist to control the model setup (Tab. 5.1), model formulation (Tab. 5.2) and the
initial conditions for all the functional groups declared in the preset layout file (see also 1.1). The
second file control the simulation environment of the standalone model, like e.g., simulation setup
(Tab. 5.4), time settings (Tab. 5.5), and environmental forcing (Tab. 5.6).

5.2. State variables initial conditions (IC)
The initial conditions of the model state variables are declared in the file BFM_General.nml within
the bfm_init_nml namelist. As a general rule, the IC of each state variable is indicated with its
code abbreviation and 0, for example Phosphate initial value will be declared as N1p0 = 1.0.
ICs for living and non-living components of the model (see Tab. 1.1) has to be treated as follow:
• Dissolved Inorganic (IO) components require the specification of the initial value;
• Living organic (LO) components require only the initial content of carbon (e.g., P1c0, Z3c0),
while the other constituents are automatically set by the model using the optimal ratios to carbon
set in the specific namelists;
• Non-living organic (NO) components require only the initial content of carbon (e.g., R1c0) and
for those groups with more than one constituent, these are set using the Redfield ratios.

5.3. Numerical integration
Three integration schemes are implemented in the STANDALONE model. They can be chosen by the
parameter method in the namelist file standalone.nml. They are all explicit and use a time step
cutting approach to avoid negative concentrations. In detail, indicating with S(cn ) the right-hand-side
source term computed with concentration cn , they are written as:
• the Euler scheme (method=1):
cn+1 = cn + S(cn )∆t

69

5. Running the STANDALONE model
name
bio_setup

kind
integer

bfm_rstctl
out_fname
out_dir
out_title

string
string
string
string

out_delta

integer

Table 5.1.: Parameters that control
BFM_General.nml)

the

name
CalcPelagic

kind
logical

CalcBenthicFlag

integer

CalcPhytoplankton(P)
CalcPelBacteria(B)

logical
logical

CalcMesoZooplankton(Z)

logical

CalcMicroZooplankton(Z)

logical

CalcPelChemistry

logical

AssignPelBenFluxesInBFMFlag logical

AssignAirPelFluxesInBFMFlag logical

ChlDynamicsFlag

integer

description
STANDALONE configuration:
1. pelagic
2. benthic (not yet available)
3. pelagic and benthic (not yet available)
Save initial state in the output file
Name of NetCDF output file
Path to the output file
Name of the experiment in NetCDF
file
Output is saved every out_delta
timesteps
model

setup

and

output

file

(bfm_nml

in

description
Switch on off the pelagic system
Switch for the benthic model (only option 1. is
currently allowed)
1. no benthic model;
2. simple benthic return;
3. benthic organisms and intermediate
complexity nutrient regeneration;
4. benthic organisms and full nutrient
regeneration (early diagenesis)
Switch for phytoplankton (P=1,2,3,4)
Switch for pelagic bacteria (B=1)
Switch for mesozooplankton (Z=1 for Z3, Z=2
for Z4)
Switch for microzooplankton (Z=1 for Z5, Z=2
for Z6)
Switch for pelagic chemistry
Switch to turn on/off the fluxes between the
pelagic and the benthic system. Fluxes are
always computed but not integrated if the flag
is false
Switch to turn on/off the fluxes between the
atmosphere and the pelagic system. Fluxes are
always computed but not integrated if the flag
is false
Switch for chlorophyll parameterization in
phytoplankton
1. constant chl:C ratio and acclimation through
optimal light (ERSEM-II). Note: the variable
for chl content P1l is used to store optimal
irradiance in W m−2
2. Chlorophyll is a phytoplankton constituent
and variable chl:C ratio. Note: P1l is in mg
chl m−3

Table 5.2.: Parameters that control model formulation and execution (Param_paramteres in
BFM_General.nml)

70

5.4. Forcing functions
name

kind

LightPeriodFlag

integer

LightLocationFlag

integer

ChlAttenFlag

integer

description
Provide the light averaging period of the
external irradiance:
1. Instantaneous irradiance
2. Daylight average with explicit photoperiod
3. Daily average
Choose the parameterization of light location
in the discrete grid:
1. Light at the top of the cell
2. Light in the middle of the cell
3. Average Light in the cell
Choose the PAR attenuation due to Chl:
1. broadband linear attenuation
2. 3-band tabulated attenuation coefficients

Table 5.3.: Parameters that control light properties in the water column (PAR_paramteres in
Pelagic_Environment.nml)
• the Runge-Kutta scheme of 2nd order (method=2):
cn+1 = cn + 12 [S(cn ) + S(cn + S(cn )∆t)] ∆t
• the leap-frog scheme (method=3):
cn+1 = cn−1 + S(cn )2∆t with the Asselin-filter: cn = cn + γ (cn−1 − cn + cn+1 )

5.4. Forcing functions
The STANDALONE model can be forced with different kind of forcing functions to provide the
physical information for the aquatic system (e.g. temperature, irradiance, etc.). Two options are
currently available: analytical forcing functions and data from an external file. All external data
formats and time interpolations are derived from the code provided by GOTM under the GPL.

5.4.1. Analytical forcing functions
This option is activated by setting forcing_method=1 in the namelist forcings_nml (Tab. 5.5).
The simple analytical forcing functions give the user a easily accessible and quickly usable basic setup of the physical forcing fields. The basic concept behind all forcing functions is the interpolation of
some (northern hemisphere) summer and winter values by a sine wave, considering them as the two
extreme values. This simple case is used for salinity forcing data (sw and ss in Table 5.6), while for
temperature, in addition to the sine wave, it is possible to superimpose a daily excursion (constant all
over the year).
Light forcing is interpolated with a sine wave from a winter minimum value to a summer maximum
value (lw and sw in Table 5.6) to obtain the actual value of average daylight for each day of the year.
These variables are as well intended as the placeholders for instantaneous light values and average
daylight, i. e. the mean value of light energy over the light period:
1
daylength

Z τdusk

light(t)dt

τdawn

The daylength is computed by the model as a function of the latitude value given in the namelist
. In case of the instantaneous light option (ltype=1) the average daylight is distributed in a sine

71

5. Running the STANDALONE model

name
nboxes

kind
integer

indepth
latitude
longitude
maxdelt
mindelt
method

real
real
real
real
real
integer

description
Number of water volumes
(boxes)
Depth of each box (m)
Latitude of each box
Longitude of each box
Maximum timestep duration (s)
Minimum timestep duration (s)
Integration method
1. Euler forward
2. Runge-Kutta 2nd order
3. Leap-frog

Table 5.4.: Main parameters for controlling time marching and integration in the STANDALONE
(standalone_nml in Standalone.nml)

name
timefmt

kind
integer

MaxN
simdays
start

integer
integer
string
“yyyy-mm-dd hh:mm:ss”
string
“yyyy-mm-dd hh:mm:ss”

stop

description
Format of the time limits and loop:
1. MaxN only. Give number of iterations,
fake start time is used.
2. start and stop dates. MaxN calculated.
3. start and MaxN. stop time calculated.
4. simdays only. Give number of simulation days, fake start time used and MaxN
calculated.
Number of iterations in time-loop
Number of simulation days
Start date
Stop date

Table 5.5.: Main parameters for time management and integration length (time_nml in
Standalone.nml))

72

5.5. Test cases
wave over the daylight period reaching thus a maximum of twice the average daylight while being
zero at the beginning and in the end of the light period. In the case of constant light (ltype=2) the
mean daylight is smoothed out over the whole day by the factor daylength/24. In the case of
rectangular light distribution (ltype=3) the light equals the average daylight value over the whole
daylight period and is zero at night time. In order to setup a constant-light experiment, it is sufficient
to put the winter value equal to the summer value and set ltype=2.

5.4.2. Boundary forcing for atmospheric CO2
The STANDALONE model (but also any other coupled configuration that does not provide a connection with the atmosphere) has the possibility to turn on the prescription of an atmospheric CO2 concentration for computing the surface fluxes. A time series can be specified as input to the atmospheric
mixing ratio by setting the file name in Carbonate_Dynamics.nml parameters listed in Tab. 2.9.
This is activated only if pelagic CO2 is turned on in the model with the macros INCLUDE_PELCO2
in the file configuration. The atmospheric value found in Carbonate_Dynamics.nml is given
in units of partial pressure (µ atm) and not in ppm. A simple annual increase of this initial value can be
turned on by giving a percentage constant rate in the variable CO2inc found in Standalone.nml
(Table 5.6, set to 0 by default, if set to 1.0 gives the typical 1% annual increase).

5.4.3. Environmental data from file
The forcing_method=2 option allows the user to input the environmental data directly from a
file. The required input data are seawater temperature (variable ETW, oC), salinity (ESW), irradiance
(in W m−2 , converted directly into the model variable EIR) and wind velocity (EWIND, m s−1 ). Any
other parameter can be included by modifying the file standalone/external_forcing.F90.
This file is a template for the input of external forcing data and must be modified by the user. The data
input file format is an ASCII file with a first column indicating the date (yyyy-mm-dd hh:mm:ss)
and one column for each parameter separated by space(s).

5.4.4. External event data
The model allows to input event data that force the state variables to certain values. This is useful to
simualte the case of a laboratory experiment, when for instance phytoplankton is inoculated at given
time interval or nutrients are restored to certain values. The syntax of the input file is the same as the
environmental data file above and the template file standalone/external_data.F90. must
be modified by the user according to the state variables that have to be modified. The subroutine is
implemented to overwrite any obtained value of the selected state variable(s) with the input value(s)
at the prescribed time(s).

5.5. Test cases
5.5.1. STANDALONE_CO2TEST: carbonate system model
The most simple test case is STANDALONE_CO2TEST that generates the run directory
run/standalone.co2test. if the BFM is configured with bfm_configure.sh -gcd
-p CO2TEST. It simulates a simple carbonate system equilibrium using the routines developed by
OCMIP (CalcCO2Method=2 in Carbonate_Dynamics.nml). Initial conditions are set for a

73

5. Running the STANDALONE model

name

kind

description

forcing_method

integer

Choice of the external forcing functions
1. analytical forcings
2. from file
3.

interactive fluxes (not yet imple-

mented)
ltype

integer

Light forcing options:
1. instantaneous light description
2. constant light over 24h
3. rectangular light distribution (on/off)
according to photoperiod
4. Constant light averaged over the photoperiod

lw

real

Sinusoidal light intensity (boreal winter) W m−2

ls

real

Sinusoidal light intensity (boreal summer) W m−2

sw

real

Sinusoidal salinity (boreal winter)

ss

real

Sinusoidal salinity (boreal summer)

tw

real

Sinusoidal temperature (boreal winter)

ts

real

degC
Sinusoidal temperature (boreal summer) degC

tde

real

Sinusoidal temperature daily excursion
degC

ww

real

Sinusoidal wind (boreal winter) m s−1

ws

real

Sinusoidal wind (boreal summer) m s−1

CO2inc

real

Linear increase in CO2 air partial pres-

forcing_file

char

sure [% per year]
Filename for external forcings (forcing_method = 2)
use_external_data
data_file

logical
char

Read external data (user defined)
Filename for external data

Table 5.6.: Parameters for the external forcing functions (forcings_nml in Standalone.nml

74

400

Atmospheric CO2 Partial Pressure[bottompoint=0]

390
380
370
360
350
340
3301980

1984

1988

1992

1996
time

2000

2004

2008

pH[oceanpoint=0]

8.30

pH[oceanpoint=0] (-)

8.25
8.20
8.15
8.10
8.05
8.001980

1984

1988

1992

1996
time

2000

2004

2008

Dissolved Inorganic Carbon[oceanpoint=0] (umol/kg)Dissolved Inorganic Carbon[oceanpoint=0] (umol/kg)

Atmospheric CO2 Partial Pressure[bottompoint=0] (uatm)

5.5. Test cases

2280

Dissolved Inorganic Carbon[oceanpoint=0]

2260
2240
2220
2200
21801980
2280

1984

1988

1992

1996
time

2000

2004

2008

Dissolved Inorganic Carbon[oceanpoint=0]

2260
2240
2220
2200
21801980

1984

1988

1992

1996
time

2000

2004

2008

Figure 5.1.: Examples from the CO2TEST simulation.
Top panels:
atmospheric CO2
partial pressure from the NOAA marine boundary layer data at 33N
(http://www.esrl.noaa.gov/gmd/ccgg/mbl/) and DIC with temperature set to the
constant value of 25 oC (default run). Bottom panels: pH and DIC values obtained with
seasonally varying temperature (tw=14, ts=27 in Standalone.nml)

75

5. Running the STANDALONE model
Mediterranean station with high alkalinity value and output is created every 60 days. The system is
forced with observed atmospheric CO2 increase over the period 1979-2009. The default simulation
prescribes a constant temperature value and it is interesting to note the differences shown in Fig. 5.1
when seasonal variability of temperature is turned on. The seasonal cycle in temperature can be prescribed by changing the winter and summer temperature minima and maxima of the analytical forcing
functions (Sec. 5.6).

5.5.2. STANDALONE_PELAGIC: the full pelagic BFM
The standard test case for the full BFM is created in directory run/standalone.pelagic and
it is the default preset when the code is generated and compiled without any option. It simulates
the seasonal cycle of a well-mixed temperate coastal sea, with mean depth 5 m. This example is
to be considered as a template for other applications, as it is not meant to describe any real-world
application. The simulation starts in January 2000 and terminates in 2010 with a maximum time step
of 8640 s in the adaptive Runge-Kutta solver and with the output stored every 30 days (see namelists
in Tab. 5.5, 5.4 and 5.1). The time evolutions of temperature, irradiance, nutrients (the system is
sustained with regenerated nutrients only, therefore nitrate plot is logarithmic and goes very close to
0), chlorophyll and gross primary production are shown in Figs. 5.2-5.4 for comparison.

5.5.3. STANDALONE_SEAICE: sea ice biogeochemistry
This configuration is the standard preset for the sea ice model described in Chap. 3. It describes the
evolution of a sea ice season using external forcing data from a sea ice model (the data are copied
from the directory inputs/seaice into the running directory by the configurator). The forcing functions
and model output are fully described in Tedesco and Vichi (2014).

5.6. Output visualization
Output values are stored in netCDF according to the CF convention (http://www.cgd.ucar.
edu/cms/eaton/cf-metadata/). The STANDALONE test cases output can be visualized with
the all-round ncview software (http://meteora.ucsd.edu/~pierce/ncview_home_
page.html), with the Python utility PyNcview https://sourceforge.net/projects/
pyncview/ or with ncbrowse. Matlab users can use the bfm_ncload.m utility available in the
tools directory. The figures presented in Sec. 5.5 have been produced with pyncview.

76

5.6. Output visualization

30

temperature[oceanpoint=0]

450

Irradiance[oceanpoint=0]

Irradiance[oceanpoint=0] (uE/m²/s)

temperature[oceanpoint=0] (C)

400

25

20

15

10

350
300
250
200
150
100
50

5

2001
time

2002

0

2001
time

2002

Figure 5.2.: Sample output of the STANDALONE_PELAGIC test case: temperature (left) and irradiance (right)

77

5. Running the STANDALONE model

Phosphate[oceanpoint=0] (mmol P/m³)

Phosphate[oceanpoint=0]
0.14
0.12
0.10
0.08
0.06
0.04
0.02
0.00

2002

2004

time

2006

2008

Nitrate[oceanpoint=0]

Nitrate (mmol N/m³)

10-1

2002

2004

time

2006

2008

Figure 5.3.: Sample output of the STANDALONE_PELAGIC test case: phosphate (top) and nitrate
(bottom)

78

0.7

Chlorophyll-a[oceanpoint=0]

0.6
0.5
0.4
0.3
0.2
0.1
0.0

Gross Primary Production (mg C/m³/d)

Chlorophyll-a[oceanpoint=0] (mg Chl/m³)

5.6. Output visualization

14

2002

2004

time

2006

2008

Gross Primary Production[oceanpoint=0]

12
10
8
6
4
2

2002

2004

time

2006

2008

Figure 5.4.: Sample output of the STANDALONE_PELAGIC test case: Chlorophyll-a (top) and gross
primary production (bottom)

79

5. Running the STANDALONE model

80

6. Structure of the code
6.1. Coding rules
The BFM is written in plain F90 and it can generally be compiled with any F90 (and above) compiler. There are no fixed coding rules apart from standard clean programming. For historical reasons,
the subroutine and symbols of the BFM biogeochemistry module are indicated with C-style capitalized long-names (e.g. EcologyDynamics), even if F90 does not distinguish capital letters. Constants
and main memory arrays are generally indicated with capital names, to highlight their global static
behavior. It is recommended to follow this style when possible.

6.2. Memory layout and variable definition
The BFM is structured similarly as the original ERSEM used to be. Pelagic variables are defined as
virtually 3-dimensional, while benthic and sea ice variables are 2-dimensional. Only the sea-points are
stored in memory, and thus the BFM has no knowledge of any land-sea mask or spatial relationships
between grid-points in general. The number of bottom points is the same as the number of surface
points, therefore the 2-D arrays can also be used for the surface variables and/or surface fluxes. These
definitions are also maintained when the model is applied to a 1-D vertical setup. In that case, the
pelagic variables have the spatial dimension of the number of vertical layers while the benthic and
surface variables have dimension 1 (the bottom and surface levels, with the index depending on the
definitions of the physical model). In a 0-D setup (as in the STANDALONE test cases described in
Chap. 5), the same grid-point is at the same time a bottom, surface and ocean layer.
Fig. 6.1 illustrates the memory layout of pelagic variables. They are gathered in a single 2dimensional array holding the number of different variables and the number of ocean grid-points:
D3STATE(NO_D3_BOX_STATES,NO_BOXES)
and the same is done for the 2-dimensional variables (defined for the benthic and sea ice ecosystem
models):
D2STATE(NO_D2_BOX_STATES,NO_BOXES_XY)
Each variable memory is accessed by means of a pointer to the corresponding row containing all
the grid-point values. This pointer array has the same naming convention of ERSEM. The rows
are numbered with named constants starting with the letters pp (ppP1c for phytoplankton carbon,
ppN1n for nitrate or ppR1p for dissolved organic phosphorus, etc.). When the BFM is coupled to
a 3-D model, then a mapping function must be taken into account as shown in Fig. 6.1 in order to
compute the transport terms.
In addition to this direct mechanism, the memory can be accessed by means of group functions.
These functions have been created to allow the addition of any number of sub-groups in the same
trophic compartment, such as diatoms, nanoflagellates, picoplankton, etc. in the group of phytoplankton. The group function is in this case called PhytoPlankton(no_of_subgroups,
no_of_components), and allows to access the memory P1c by writing

81

6. Structure of the code

Figure 6.1.: Layout of the main memory array for the pelagic variables (D3STATE) and schematic of
the remapping into a three-dimensional ocean grid for transport processes.
PhytoPlankton(iiP1,iiC)
The model state variables and components are defined in the model layout file of each preset (Sec.
4.2) which will be detailed in Chap. 7. It is now only sufficient to list the constants indicating the constituents, namely iiC, iiN, iiP and iiSi. For each defined group there is a named constant starting
with the letters ii which holds the total number of existing subgroups (e.g. iiPhytoPlankton).
A typical example of the use of groups is to compute the total amount of food available to microzooplankton (from MicroZoo.F90):
do i = 1 ,iiPhytoPlankton
PPYc(:,i) = p_paPPY(zoo,i)*PhytoPlankton(i,iiC)* &
MM_vector(PhytoPlankton(i,iiC), p_minfood(zoo))
rumc = rumc + PPYc(:,i)
rumn = rumn + PPYc(:,i)*qncPPY(i,:)
rump = rump + PPYc(:,i)*qpcPPY(i,:)
end do
Compare this code with the equation for food availability and capture efficiency presented in Sec.
2.4.1 of the BFM description and it will show how it is simple to cycle over all the phytoplankton
variables, independently of the number and type of the defined subgroups.
The other basic arrays of the memory are the ones containing the time rates of change due to each
physiological or ecological interaction. The rate of mass exchange between two state variables is an

82

6.3. The BFM flow-chart
element of an array D3SOURCE, if it is a production term, and D3SINK if it is a loss term, both
with dimensions NO_D3_BOX_STATES × NO_D3_BOX_STATES × NO_BOXES (Fig. 6.2). It is
clear that the two matrices are symmetrical and a predation (negative) rate from variable P1c to Z4c
is equivalent to an ingestion rate (positive) of Z4c over P1c. All the rates are defined positive by
convention, so that to obtain the net rate of change for variable of row j, it is sufficient to compute
sum(D3SOURCE(j,:,:)-D3SINK(j,:,:))
Fluxes which involve a variable not defined in the model, such as oxygen production from water
during photosynthesis can be taken into account by using the diagonal of the matrix. There is a
positive flow from O2o to itself computed stoichiometrically from the gross carbon production rate
and inserted in (D3SOURCE(O2o,O2o,:), cf. Fig. 6.2).
The same definitions apply for the 2-dimensional system, and the arrays are called D2SOURCE and
D2SINK with dimensions NO_D2_BOX_STATES × NO_D2_BOX_STATES × NO_BOXES_XY.
This separation between positive and negative terms allow to use more sophisticated positivedefinite integration schemes, but it may result in an excessive memory usage. Thus, the default setting
of BFM make use of a shrank version of D3SOURCE (NO_D3_BOX_STATES × NO_BOXES), which
is the result of the addition of all the terms after the combination of the two symmetrical matrices into
one that contains all the source (upper-left part) and sink (lower-right part) terms. It is possible to use
the two separated matrices for source and sinks terms through the inclusion of the compilation macro
EXPLICIT_SINK in the preset configuration file.
No direct mass exchanges are allowed between variables from different realms (like for instance
between pelagic and benthic variables), and therefore those fluxes are defined as diagonal fluxes
each one within its own realm taking into account any change of units. Specialized functions
have been written to take into account all the possible kind of exchanges and are collected in
General/FluxFunctions.F90.

6.3. The BFM flow-chart
The BFM core will be described in its STANDALONE configuration. The main program (Fig. 6.3)
consists of 3 subroutines, one for initialization, one for the time-loop and one for closing all active
processes. The same structure is likely to be maintained also for all the coupled configurations, but
the calls are going to be distributed differently according to the structure of the general circulation
model.

6.3.1. Initialization procedures
The initialization of the BFM requires 4 different phases:
1. Initialization of geometric domain (number of boxes/gridpoints).
This part depends on the configuration:
in the STANDALONE case, this is done in
standalone/standalone.F90[init_standalone] while in the other coupled
configuration there are specific interface subroutines;
2. Initialization of the set-up via namelist and definition of variable information
(share/api_bfm.F90[init_bfm] for the STANDALONE case);
3. Definition of variable information and allocation/initialization of variables;
4. Initialization of the NetCDF output.

83

6. Structure of the code

Figure 6.2.: Schematic of the D3SOURCE array collecting the source terms.

6.3.2. Time marching
The time marching is driven by the host hydrodynamical model, or in the case of the STANDALONE
model, by a specific built-in routine. The phases of the time stepping shown in Fig. 6.3 are the
following:
1. update of the physical forcing functions needed by the biogeochemical model (temperature, irradiance, wind, etc.). This is done by subroutine envforcing_bfm, which in the case of the
STANDALONE model computes the analytical values based on the time of the day/year. When
coupled with an ocean model, this routine is the contact point between the two models. Note
that information may flow in two directions because the light attenuation coefficient might also
be passed back to the physical model via this routine.
2. computation of the reaction term for the BFM module. This is done in EcologyDynamics,
which is further explained in the next section.
3. numerical integration of the biogeochemical terms. In the Standalone model there are three
different numerical schemes that can be used (see Sec. 5.3). Note that EcologyDynamics
can be called more than one time by this routine when the integration scheme is of higher-order.
4. preparation and writing of the NetCDF output, routines calcmean_bfm and save_bfm (the
names of the routines can be different for each coupling according to the ocean model structure
and naming conventions).

84

6.3. The BFM flow-chart

Figure 6.3.: Invocation tree of the Standalone main.

85

6. Structure of the code

Figure 6.4.: Invocation tree for PelagicSystemDynamics

86

6.3. The BFM flow-chart

6.3.3. Computation of pelagic reaction terms
The pelagic components are computed in PelagicSystemDynamics, which invokes the routines shown in Fig. 6.4. It first sets the diagnostic values for total chlorophyll (sum of all the
chlorophyll components in phytoplankton) and for the nutrient:carbon ratios in all the components
(PelGlobalDynamics) and for oxygen saturation and surface fluxes with the atmosphere.
The routines for the computation of light availability are only used when the parameterization of
optimal irradiance is activated (as in ERSEM-II, Ebenhöh et al., 1997), see Sec. 2.2.7, which is done
in Param_parameters (Tab. 5.2) by setting ChlDynamicsFlag=1.
The standard groups, phytoplankton, meso- and micro-zooplankton and bacteria are computed
in PhytoDynamics, Meso- MicroZooDynamics and PelBacDynamics, respectively. The
same routine is modularly used for each sub-group, as
!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
! Compute phytoplankton dynamics
!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
do i =1,iiPhytoPlankton
if ( CalcPhytoPlankton(i)) call PhytoDynamics( i )
end do
...
!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
! Compute MesoZooPlankton dynamics
!-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
do i =1,iiMesoZooPlankton
if ( CalcMesoZooPlankton(i)) call MesoZooDynamics( i )
end do

PhytoDynamics has the declaration structure depicted in Fig. 6.5. This is an example
of how the memory layout presented in Sec. 6.2 is accessed by each submodel of the BFM.
To refer to the pointer of each phytoplankton component, one can use the placeholder function
ppPhytoPlankton(i,j), where i is the phytoplankton subgroup number (up to 4 in the reference
model) and j is the component (c, n, p, etc). Using ppPhytoPlankton(1,iiC) is equivalent to
writing ppP1c.
The function flux_vector is invoked many times by PhytoDynamics (Fig. 6.5) because it is
used to fill in the source/sink term arrays described in Sec. 6.2. This function is used throughout the
BFM modules. In case of a mass exchange term such as the excretion of dissolved carbon we have
call flux_vector(iiPel,ppphytoc,ppR1c,rr1c)
where the flux is pelagic (iiPel), from phytoplankton carbon to dissolved organic carbon
(ppphytoc, ppR1c) and the value is rr1c, which will be assigned both to D3SINK and
D3SOURCE in the appropriate element. The oxygen production rate is instead assigned as
call flux_vector(iiPel,ppO2o,ppO2o,rugc/MW_C)
which concerns a pelagic flux (iiPel), on the oxygen diagonal (from ppO2o to ppO2o) and MW_C
is the mass weight of carbon to convert gross primary production (rugc). Since the last term is
positive, this flux will be assigned to D3SOURCE.
The same methodology is applied for all the other pelagic fluxes and also in the benthic and sea
ice system. After the computation of all the reaction terms for pelagic organisms, there is the call

87

6. Structure of the code

Figure 6.5.: Declaration scheme of PhytoDynamics.
to pelagic chemistry (PelChemDynamics, see Fig. 6.4). This part computes all the parameterized bacteria-mediated biochemical processes such as nitrification and denitrification, re-oxidation of
reduction equivalents and dissolution of biogenic silica.

88

7. Design model layout components
7.1. Layout syntax and modular structure
The BFM is conceived in a modular way and the number of state variables and biogeochemical basic components can be modified according to the modelling needs. In fact, each
state variable is not defined directly in the code, but through a model template layout
found in the PRESET directories (Sec. 4.2, $BFMDIR/build/configurations). All the
ancillary F90 memory layout files (AllocateMem.F90, INCLUDE.h, ModuleMem.f90,
set_var_info_bfm.f90) are automatically generated from this meta-data structure by way of
a perl script ($BFMDIR/build/script). The prototypes of these files are found in directory
$BFMDIR/BFM/proto and filled in with the user-defined state variables and diagnostics.
The structure of the model template given in Tab. 7.1 is divided in sections with specific headings.
Sections can be repeated many times in the file. For clarity, it is suggested to follow a conceptual
construction of the model, by grouping all the state variables and diagnostics as done in the standard
released structure.
3d-state: this section is used for all the state variables, which are defined as being virtually threedimensional. This means that they exist in the whole domain, be it 0-, 1- or 3-dimensional. An
example of state variable definition is shown in Box 7.1.
According to the theoretical model description, variables are defined as functional families and
identified by their main biochemical components, which are indicated according to the classical ERSEM syntax (Blackford and Radford, 1995). State variable definitions require 3 pieces of
information, separated by colons: code symbol with components, long name and units. Components are indicated by one single letter and currently the following are reserved: c=carbon,
n=nitrogen, p=phosphorus, s=silica, l=chlorophyll, f=iron. Any other new component (either atomic or molecular) can be added, provided that the user include the relative dynamics
in an appropriate piece of code. The syntax for a multi-component variable is detailed below varname[cnp...], and it is also possible to define group of state variables, such as
PhytoPlankton or PelDetritus as shown in Box 7.1. This feature gives the advantage to
cycle over all the group components with one single index (see the example in Sec. 6.2).
2d-state: all the features described above are also valid for the state variables belonging to the
benthic or sea ice 2-D system.
1d-variable: this section defines the variable without a spatial component (more appropriately defined as 0-D variables). These are only function of time and defined as scalars.
2d-variable: this section includes the surface variables, originally identified as benthic variables.
However, this section is also used to define surface variables such as wind velocity. In a 1D
coupled model, this is a scalar, while it is a one-dimensional array holding all the surface (bottom)
points in a coupled 3D model.

89

7. Design model layout components

...
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
PELAGIC STATE VARIABLES (volume)
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=3d-state
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
State Variable(s) for Nutrients
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=N1p :
Phosphate : mmol P/m3
N3n :
Nitrate : mmol N/m3
N4n :
Ammonium : mmol N/m3

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
State Variable(s) for Phytoplankton
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=group PhytoPlankton[cnpls] (PPY) : mg C/m3 : mmol N/m3 : mmol P/m3 :
mg Chl/m3 : mmol Si/m3
P1
: Diatoms
P2[-s]
: Flagellates
P3[-s]
: PicoPhytoplankton
P4[-s]
: Large Phytoplankton
end
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
State Variable(s) for Mesozooplankton
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=group MesoZooPlankton[cnp] (MEZ) : mg C/m3 : mmol N/m3 : mmol P/m3
Z3
: Carnivorous Mesozooplankton
Z4
: Omnivorous Mesozooplankton
end
end

Table 7.1.: Example of the model template file layout with the definition of pelagic state variables.

90

7.1. Layout syntax and modular structure
3d-variable: this section is used for all the ancillary diagnostic variables defined in the whole domain. Typically, this section is used for temperature, salinity, and other environmental or diagnostic variables. Environmental variables obtained from external data files or hydrodynamical
models are indicated with capital symbols starting with E (e.g. ETW = environmental temperature
of water).
The syntax of the allowed layout elements is the following
[dim]d-state [-if-exist]
GROUP
VARNAME
end
[dim]d-variable [-if-exist]
VARNAME
ARRAYNAME
QUOTUM
end
[dim]d-flux [-if-exist]
FLUXNAME
end

• GROUP, like PhytoPlankton
group name[constituents] (acronym) : constituent_units
VARNAME
end

• VARNAME, like N1p or P2[-s] (if present, the constituent within brackets are removed from
the group list of elements)
name[-constituents] : description : units

• ARRAYNAME, a variable that is defined for each element of a group, like the specific uptake
for each phytoplankton suPPY(PhytoPlankton)
name(groupname) : description : units

• QUOTUM, like qncPPY(PhytoPlankton) (quotum N:C for each phytoplankton group)
q+constituent+constituent+acronym() : description : units

• FLUXNAME, the named entity for a user-defined flux like ruPPYc (Sec. 8.3)
name = FLUX : description : units

• FLUX, an expression combining several variables to store fluxes between components (Sec. 8.3)
(var[+|-]var)[<-|->](var[+|-]var)

with the following fundamental elements

91

7. Design model layout components
• name: an alphanumeric string starting with a non-numeric character
• constituents: one or several from the list: c (Carbon), n (nitrogen), p (phosphorus), s (silicon), f (iron), l (chlorophyll). They have to be written in lower case and without any separation
• description: phrase describing the variable, quota or flux
• units: units in which the variable, quota or flux are measured
• constituents_units: units separated by “:”. Have to be same number of units as constituents have been specified
• groupname: name of a group, also name of the function to access the group components
• var: variable name (e.g. R6c) or first letter of group with constituents for naming all (e.g. P.c)
• acronym: acronym representing the group. Currently accepted acronyms are: PPY
(PhytoPlankton), PBA (PelBacteria), MEZ (MesoZooPlankton), MIZ
(MicroZooPlankton), OMT (PelDetritus), CAR (Inorganic)
• dim: dimension of the layout elements inside. Allowed values are: 1, 2 and 3
The configuration can be also modified at compilation time by means of the set of macros listed in
Sec. 4.2. The statement -if-exist found after one of the layout declarations indicates that the
following set of variables will be defined only when the macro is included in the compilation:
3d-variable -if-exist INCLUDE_PELCO2
DIC
: Dissolved Inorganic Carbon : umol/kg
CO2
: CO2(aq) : umol/kg
pCO2
: Oceanic pCO2 : uatm
HCO3
: Bicarbonate : umol/kg
CO3
: Carbonate : umol/kg
Ac
: Alkalinity : umol eq/kg
pH
: pH : end

7.2. Example 1. Adding a subgroup
This is the most simple example, because it involves a small number of changes. Let us assume we
want to add another mesozooplankton component to the zooplankton group , such as macrozooplankton (as done for instance by Berline et al., 2011). It is also assumed that this group will have exactly
the same dynamics as the other two subgroups, but different parameter values. It is first suggested to
create a new preset by copying the directory STANDALONE_PELAGIC to a new name
cd $BFMDIR/build/configurations
cp -r STANDALONE_PELAGIC STANDALONE_PELAGIC_Z

and then edit the appropriate line in the template file layout by adding the new variable name and
long name:

92

7.2. Example 1. Adding a subgroup
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
State Variable(s) for Mesozooplankton
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=group MesoZooPlankton[cnp] (MEZ) : mg C/m3 : mmol N/m3 : mmol P/m3
Z3
: Carnivorous Mesozooplankton
Z4
: Omnivorous Mesozooplankton
Z2
: Macrozooplankton
end

It is also suggested to edit the name of the run directory where the namelists will be copied. This is
done in the configuration file $BFMDIR/build/configurations/configuration and change
EXP

= ’standalone.macrozoo’

Now the configuration can be generated by compiling the new preset
cd $BFMDIR/build
./bfm_configure.sh -gcd -p STANDALONE_PELAGIC_Z

and the script will take care of the automatic generation of the new memory layout with the required
named constants, pointers, initial value, run directory and namelists. A new Z2 element will be added
to the group MesoZooPlankton (in this case Z2 will be the third element as it comes after Z3 and
Z4). Particularly, a new column will be added in the mesozooplankton namelist and filled in with
zero. A warning is issued for each added column or row of parameters:
WARNING:
WARNING:
WARNING:
...
WARNING:
WARNING:

(2 < 3) adding zero values to element p_q10 in namelist MesoZoo
(2 < 3) adding zero values to element p_srs in namelist MesoZoo
(2 < 3) adding zero values to element p_sum in namelist MesoZoo
(2 < 3) adding new predator term p_paMIZ(3,:) in namelist MesoZoo
(2 < 3) adding new predator term p_paPPY(3,:) in namelist MesoZoo

Check the updated file tmp/STANDALONE_PELAGIC_Z/ModuleMem.F90 and search for the
new variable Z2 to see the added pieces of code.
The final steps before running the new model are the addition of the proper parameter values
in place of zeros in the new column/rows corresponding to Z2 of MesoZoo_parameters in
run/STANDALONE_PELAGIC_Z/Pelagic_Ecology.nml
&MesoZoo_parameters
!
p_q10
p_srs
p_sum
p_vum
p_puI
p_peI
p_sdo
p_sd
p_sds
p_qpcMEZ
p_qncMEZ
p_clO2o
!
!
Z3
p_paPPY(1,:)

=
=
=
=
=
=
=
=
=
=
=
=

=

Z3
2.0
0.01
2.0
0.025
0.6
0.3
0.01
0.02
2.0
1.67d-3
0.015
30.0
P1

Z4
2.0
0.02
2.0
0.025
0.6
0.35
0.01
0.02
2.0
1.67d-3
0.015
30.0
P2

Z2
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
P3

P4

0.0

0.0

0.0

0.0

93

7. Design model layout components
!
!
!
!
!
!
!
!
!
!
!
!
!

Z4
p_paPPY(2,:)
Z3
p_paMIZ(1,:)
Z4
p_paMIZ(2,:)
Z3
p_paMEZ(1,:)
Z4
p_paMEZ(2,:)
Z2
p_paMIZ(3,:)

=

1.0
Z5

0.0
Z6

0.0

=

0.0

0.0

=

1.0
Z3

0.0
Z4

Z2

=

1.0

1.0

0.0

=

0.0
Z5

1.0
Z6

0.0

=

0.0
P1

0.0
P2

P3

P4
0.0

Z2
p_paPPY(3,:)

=

0.0
Z3

0.0
Z4

0.0
Z2

Z2
p_paMEZ(3,:)

=

0.0

0.0

0.0

0.7

/

and to turn on the flag CalcMesoZooPlankton(3)=.TRUE. in Param_parameters and add
initial value(s) by adding a line Z2c0= in bfm_init_nml, both found in file
run/STANDALONE_PELAGIC_Z/BFM_General.nml.

7.3. Example 2. Adding a biochemical component and a group
The BFM application in the global ocean is called PELAGOS (PELAgic biogeochemistry for Global
Ocean Simulations) and the climatological simulations were described in Vichi et al. (2007a). A major requirement to simulate the global ocean ecosystem is the introduction of iron control on phytoplankton growth (Sec. 2.2.8), particularly to capture the high-nutrient low-chlorophyll regions. The
example below shows how the iron component of phytoplankton and detritus have been added to the
model structure template. Additional coding from the user is then required to parameterize the fluxes
between these components and the regulations of plankton physiology (this is already included in
the GYRE_BFM preset that requires the usage of the NEMO model and activated with the keyword
INCLUDE_PELFE).
To add an iron constituent to a group, it is sufficient to specify a unique single-character constant
between the square brackets in the layout file as shown below:
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#
State Variable(s) for Phytoplankton
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
group PhytoPlankton[cnplsf] (PPY) : mg C/m3 : mmol N/m3 : mmol P/m3 :
mg Chl/m3 : mmol Si/m3 : umol Fe/m3
P1
: Diatoms
P2[-s] : Flagellates
P3[-s] : PicoPhytoPlankton
end

94

7.4. Example 3. Removing components from zooplankton
It is also necessary to provide units separated by colons, in this case µ mol Fe m−3 . Remember
that the [-s] syntax in some phytoplankton groups indicates that this constituent, although being a
component of the group vector, is not computed for that specific sub-group.
The preprocessing script will take care during compilation to generate an adequate entry in the main
memory module and the definitions of the named constants (ppP1f, . . . ) and pointers (P1f, . . . ) to
access the new variables from the model code.
New groups can also be defined in order to use the facilities shown in Sec. 6.3.3. All the detritus
variables can be grouped together with the following syntax:
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#
State Variable(s) for Detritus (Biogenic Organic Material)
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=group PelDetritus[cnpsf] (OMT) : mg C/m3 : mmol N/m3: mmol P/m3:
mmol Si/m3 : umol Fe/m3
R1[-s]
: Labile Dissolved Organic Matter
R2[-npsf] : Semi-labile Dissolved Organic Carbon
R3[-npsf] : Semi-refractory Dissolved Organic Carbon
R6
: Particulate Organic Matter
end

7.4. Example 3. Removing components from zooplankton
For certain purposes it is not needed to carry out the computation of variable nutrient:carbon quota
in selected Functional Groups (for instance, in the case of the new macrozooplankton group added in
example 1). This can be done to save computational resources in the case of low-resolution long-term
global ocean simulations or with high-resolution coastal models. A possible example is the assumption of a fixed quota in zooplankton which allows one to reduce the number of standard components
of up to 8 state variables. The zooplankton group definitions are equal to the one shown in the code
of Box 7.1 for MesoZooplankton but with the removal of the n and p constituents:
#−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−
#
S t a t e V a r i a b l e ( s ) f o r Mesozooplankton
#−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−=−
g r o u p M esoZ ooP lankton [ cnp ] (MEZ) : mg C / m3 : mmol N / m3 : mmol P / m3
Z3[−np ]
: C a r n i v o r o u s Mesozooplankton
Z4[−np ]
: Omnivorous M e s o z o o p l a n k t o n
end
The additional requirement is the change of value in parameter check_fixed_quota from 0 to 1
in the Param_parameters namelist. Test experiments have shown that the simulation times can
be reduced up to 30%. The physiological assumption is that zooplankton release in inorganic form
the nutrients in excess of the given fixed quota.

95

7. Design model layout components

96

8. Model outputs
8.1. Output
The variables that can be stored in the output file are the ones defined in file layout. All the
state variables can be stored (provided that they have a long name and units), both as instantaneous
values and as means over the chosen output period. This behavior is controlled via the namelist
bfm_save_nml contained in the BFM_General.nml file when running the STANDALONE
model. Other couplings may have different names. Additionally, a set of other variables can be saved,
such as the ones obtained from the physical model and other diagnostics such as total chlorophyll,
nutrient quota or aggregated rates.

8.2. Restart
The length of each BFM simulation is controlled by the physical transport model except for the STANDALONE model which has its own time marching. It is common practice in modeling to split the
simulation in several chunks and to continue the simulation from the end values of each state variable
at the end of a previous run. The BFM stores this information in a NetCDF file which is always generated at the end of each simulation. There is no automatic control on the date (neither it is stored as
an attribute inside the NetCDF file), so it is totally up to the user to check that the initial restart file is
the proper one. The input/output restart file names and the restart switch (bfm_init) are set in the
BFM_General.nml file within the namelist bfm_init_nml.

8.3. Aggregated diagnostics and rates
All the ancillary 2D and 3D variables (like wind velocity or temperature, for instance) and diagnostics
are grouped in two memory arrays, D3DIAGNOS and D2DIAGNOS, which are accessed by means of
pointers. The dimensions of these arrays change according to the number of model elements in the
layout file, and any diagnostic rate can be easily defined thanks to the automatic generation of the
memory layout. Some fluxes are already defined by default, such as the gross/net primary production
rates, total respiration rates and (de)nitrification rates.
Any exchange of mass between every state variable can be tracked by means of a special
syntax that allows the user to access and combine the elements of the D3(D2)SOURCE and
D3(D2)SINK arrays where all the rates are stored. Some examples are already included at the
end of the STANDALONE_PELAGIC/layout file and are activated with the compilation macros
INCLUDE_DIAG:
#
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=#

DIAGNOSTIC RATES

97

8. Model outputs
#
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=3d-flux -if-exist INCLUDE_DIAG INCLUDE_PELCO2
ruPPYc=P.c <- O3c : Gross Primary Production : mg C/m3/d
resPPYc=P.c -> O3c : Phytoplankton Respiration
: mg C/m3/d
resPBAc=B.c -> O3c : Bacterial Respiration
: mg C/m3/d
resZOOc=Z.c -> O3c : Zooplankton Respiration
: mg C/m3/d
end
3d-flux -if-exist INCLUDE_DIAG
ruPPYn=P.n <- N3n+N4n: Net Nitrogen Uptake:mmol N/m3/d
ruPPYp=P.p <- N1p : Net Phosphate Uptake : mmol P/m3/d
ruPPYs=P.s <- N5s : Net Silicate Uptake : mmol Si/m3/d
exPPYc=(P.c->R1c+R2c+R3c+R6c): C excretion from phytoplankton :mg C/m3/
d
ruZOOc=(Z.c<-P.c+B.c+Z.c): Gross Secondary Production:mg C/m3/d
resPELo=(O2o->*): Pelagic Oxygen Demand : mmol O2/m3/d
remZOOn=(N4n<-B.n+Z.n): Zooplankton N Remineralization : mmol N/m3/d
remZOOp=(N1p<-B.p+Z.p): Zooplankton P Remineralization : mmol P/m3/d
remPBAn=(N4n<-B.n): Bacterial N Uptake/Remineralization : mmol N/m3/d
remPBAp=(N1p<-B.p): bacterial P Uptake/Remineralization : mmol P/m3/d
end

The names used for the diagnostics are derived from the short names of the groups defined above
in the layout file (see Sec. 7), preceded by the letter r for rate and the identifier of the chemical
constituent. This is not mandatory, but it is suggested to follow the convention derived from the
original ERSEM code (Blackford and Radford, 1995) to easily identify the type of rate. The first
line, for instance, defines the combined rate of gross primary production ruPPYc (rate of uptake of
pelagic phytoplankton PPY expressed in carbon), which is the sum of all the carbon fluxes to all the
phytoplankton groups (notice the syntax using the wildcard “.” to indicate all sub-groups).
An example of aggregated flux is the carbon excretion from phytoplankton exPPYc or gross secondary production ruZOOc, which is defined as the sum of carbon ingestion by all zooplankton
groups. The diagnostics also define the long_name and units attributes for the NetCDF output.
Diagnostics can be included in the output just by adding the name to the bfm_save_nml namelist
in BFM_General.nml.
Boundary fluxes are defined by default for each pelagic state variable. The variables
PELSURFACE(NO_D3_BOX_STATES,NO_BOXES_XY)
PELBOTTOM(NO_D3_BOX_STATES,NO_BOXES_XY)

are also pointers to the main diagnostic memory D2DIAGNOS and store the surface and bottom
boundary flux for any state variable. They can be accessed by means of specific pointers allocated
in the generated subroutine BFM/General/AllocateMem.F90 (notice that the integer constants
change at every code generation and are only reported here as an example of the code):
PELSURFACE
jsurO2o =>
jsurN1p =>
jsurN3n =>
jsurN4n =>

=> D2DIAGNOS(17+1:17+50,:); PELSURFACE=ZERO
D2DIAGNOS(17+ppO2o,:); jsurO2o=ZERO
D2DIAGNOS(17+ppN1p,:); jsurN1p=ZERO
D2DIAGNOS(17+ppN3n,:); jsurN3n=ZERO
D2DIAGNOS(17+ppN4n,:); jsurN4n=ZERO

The boundary rates (always in units per day) are accessed in the code by means of both the direct

98

8.4. Mass conservation
value PELSURFACE(ppO2o,:) or the pointer jsurO2o(:), but can be added to the output in the
bfm_save_nml section of BFM_General.nml only with the pointer name (e.g. jsurO2o).

8.4. Mass conservation
The total mass in the system (mass units m−2 ) can be checked by means of built-in diagnostic variables. There is a variable for each biogeochemical component (C, N, P, Si, etc.) both for the pelagic
and the benthic systems. To activate the diagnostic it is necessary to add the following lines to the
output namelist in the var_save section:
var_save = ’totpeln’,
’totpelp’,
’totpels’

The actual computation is done in General/CheckMassConservation.F90. Note that for
technical reasons the total mass variables are multidimensional but the value is stored in the first
element of the array. When checking a coupling with an OGCM, it is also important to turn off any
additional boundary flux (e.g. dilution, atmospheric deposition, river input, etc.) because they are not
taken into account by the standard BFM.

99

8. Model outputs

100

Bibliography
Andersen, T., Elser, J. J., Hessen, D. O., 2004. Berline, L., Stemmann, L., Vichi, M., Lombard,
Stoichiometry and population dynamics. Ecol.
F., Gorsky, G., 06 2011. Impact of appendicLett. 7, 884–900.
ularians on detritus and export fluxes: a model
approach at dyfamed site. J. Plankt. Res. 33 (6),
Anderson, T. R., 2005. Plankton functional type
855–872.
modelling: running before we can walk? J.
Plankt. Res. 27, 1073–1081.
Bidle, K. D., Azam, F., 2001. Bacterial control of silicon regeneration from diatom detriArcher, D. E., Johnson, K. S., 2000. A Model of
tus: Significance of bacterial ectohydrolases
the iron cycle in the ocean. Glob. Biogeochem.
and species identity. Limnol. Oceanogr. 46 (7),
Cy. 14, 269–279.
1606–1623.
Arrigo, K. R., Dieckmann, G., Gosselin, M.,
Robinson, D. H., Fritsen, C. H., Sullivan, Blackford, J., Radford, P., 1995. A structure and
methodology for marine ecosystem modelling.
C. W., 1995. High resolution study of the
J. Sea Res. 33 (3-4), 247–260.
platelet ice ecosystem in McMurdo Sound,
Antarctica: biomass, nutrient, and production
Blackford, J. C., Allen, J. I., Gilbert, F. J., 2004.
profiles within a dense microalgal bloom. Mar.
Ecosystem dynamics at six contrasting sites: a
Ecol. Prog. Ser. 127, 255–268.
generic modelling study. J. Mar. Sys. 52, 191–
215.
Aumont, O., Maier-Reimer, E., Monfray, P.,
Blain, S., 2003. An ecosystem model of the
Blackford, J. C., Burkill, P. H., 2002. Planktonic
global ocean including Fe, Si, P co-limitations.
community structure and carbon cycling in the
Glob. Biogeochem. Cy. 17 (2), 1060.
Arabian Sea as a result of monsoonal forcing:
Baretta, J., Ebenhöh, W., Ruardij, P., 1995. The
the application of a generic model. J. Mar. Sys.
European Regional Seas Ecosystem Model, a
36 (3), 239–267.
complex marine ecosystem model. J. Sea Res.
Broekhuizen, N., Heath, M., Hay, S., Gurney,
33 (3-4), 233–246.
W., 1995. Modelling the dynamics of the North
Baretta-Bekker, J., Baretta, J., Ebenhoeh, W.,
Sea’s mesozooplankton. J. Sea Res. 33 (3-4),
1997. Microbial dynamics in the marine
381–406.
ecosystem model ERSEM II with decoupled
carbon assimilation and nutrient uptake. J. Sea Contois, D., 1959. Kinetics of bacterial growth:
relationship between population density and
Res. 38 (3/4), 195–212.
specific growth rate of continuous cultures.
Baretta-Bekker, J., Baretta, J., Rasmussen, E.,
Journal of general microbiology 21 (1), 40–50.
1995. The microbial food web in the European
Regional Seas Ecosystem Model. J. Sea Res. Cota, G. F., Legendre, L., Gosselin, M., Ingram,
33 (3-4), 363–379.
R. G., 1991. Ecology of bottom ice algae: I
Environmental controls and variability. J. Mar.
Behrenfeld, M. J., Prasil, O., Babin, M., Bruyant,
Syst. 2, 257–277.
F., 2004. In search of a physiological basis for
covariations in light-limited and light-saturated
photosynthesis. J. Phycol. 40, 4–25.

101

Bibliography
Dieckmann, G. S., Lange, M. A., Ackley, S. F., Geider, R., MacIntyre, H., Kana, T., 1997. A dyJennings, J. C., 1991. The Nutrient Status in
namic model of phytoplankton growth and acSea Ice of the Weddell Sea During Winter: Efclimation: responses of the balanced growth
fects of Sea Ice Texture and Algae. Polar Biolrate and chlorophyll a:carbon ratio to light, nuogy 11, 449–456.
trient limitation and temperature. Mar. Ecol.
Prog. Ser. 148, 187–200.
Doney, S. C., Lindsay, K., Caldeira, K., Campin,
J. M., Drange, H., Dutay, J. C., Follows, M., Geider, R., MacIntyre, H., Kana, T., 1998. A dyGao, Y., Gnanadesikan, A., Gruber, N., Ishida,
namic regulatory model of phytoplanktonic acA., Joos, F., Madec, G., Maier-reimer, E.,
climation to light, nutrients, and temperature.
Marshall, J. C., Matear, R. J., Monfray, P.,
Limnol. Oceanogr. 43 (3), 679–694.
Mouchet, A., Najjar, R., Orr, J. C., Plattner,
G. K., Sarmiento, J., Schlitzer, R., Slater, R., Gentleman, W., Leising, A., Frost, B., Strom, S.,
Murray, J., 2003. Functional responses for zooTotterdell, I. J., Weirig, M. F., Yamanaka, Y.,
plankton feeding on multiple resources: a reYool, A., 2004. Evaluating global ocean carbon
view of assumptions and biological dynamics.
models: The importance of realistic physics.
Deep-Sea Res. Pt. II 50, 2847–2875.
Glob. Biogeochem. Cy. 18, 3017.
Droop, M., 1973. Some thoughts on nutrient lim- Gibson, G. A., Musgrave, D. L., Hinckley, S.,
2005. Non-linear dynamics of a pelagic ecosysitation in algae. J. Phycol. 9, 264–272.
tem model with multiple predator and prey
Ebenhöh, W., Baretta-Bekker, J., Baretta, J.,
types. J. Plankt. Res. 27, 427–447.
1997. The primary production module in the
marine ecosystem model ERSEM II with em- Goldman, J., McCarthy, J., 1978. Steady
phasis on the light forcing. J. Sea Res. 38, 173–
state growth and ammonium uptake of a
193.
fast-growing marine diatom. Limnology and
oceanography 23 (4), 695–703.
Flynn, K. J., 2 2003. Modelling multi-nutrient interactions in phytoplankton; balancing simplic- Gosselin, M., Legendre, L., Demers, S., Ingram,
ity and realism. Prog. Oceanogr. 56 (2), 249–
R. G., 1985. Responses of Sea-Ice Microalgae
279.
to Climatic and Fortnightly Tidal Energy Inputs
(Manitounuk Sound, Hudson Bay). Canadian
Flynn, K. J., Marshall, H., Geider, R. J., 2001.
Journal of Fisheries and Aquatic Sciences 42,
A comparison of two N-irradiance interac5, 999–1006, doi:10.1139/CJFAS–42–5–999.
tion models of phytoplankton growth. Limnol.
Oceanogr. 46, 1794–1802.
Granskog, M. A., Kaartokallio, H., Kuosa, H.,
Thomas, D. N., Vainio, J., 2006. Sea ice in
Follows, M. J., Ito, T., Dutkiewicz, S., 2006. On
the baltic sea: A review. Estuarine, Coastal and
the solution of the carbonate chemistry system
Shelf Science 70, 145–160.
in ocean biogeochemistry models. Ocean Modelling 12 (3-4), 290–301.
Heissenberger, A., Leppard, G. G., Herndl, G. J.,
Frost, P. C., Xenopoulos, M. A., Larson, J. H.,
1996. Relationship between the intracellular in2004. The stoichiometry of dissolved organic
tegrity and the morphology of the capsular encarbon, nitrogen, and phosphorus release by a
velope in attached and free-living marine bacplanktonic grazer, Daphnia. Limnol. Oceanogr.
teria. Applied and environmental microbiology
49, 1802–1808.
62 (12), 4521–4528.
Geider, R., MacIntyre, H., Kana, T., 1996. A
dynamic model of photoadaptation in phytoplankton. Limnol. Oceanogr. 41 (1), 1–15.

102

Bibliography
Jassby, A., Platt, T., 1976. Mathematical formula- Morel, F., 1987. Kinetics of uptake and growth in
tion of the relationship between photosynthesis
phytoplankton. J. Phycol. 23, 137–150.
and light for phytoplankton. Limnol. Oceanogr.
Olsen, A., Wanninkhof, R., Trinanes, J. A., Jo21, 540–547.
hannessen, T., 2005. The effect of wind speed
Johnson, K. S., Gordon, R. M., Coale, K. H.,
products and wind speed-gas exchange rela1997. What controls dissolved iron concentrationships on interannual variability of the airtions in the world ocean? Mar. Chem. 57, 137–
sea CO2 gas transfer velocity. Tellus B 57, 95–
161.
106.
Lazzari, P., Solidoro, C., Ibello, V., Salon, S., Parekh, P., Follows, M. J., Boyle, E., 2004. ModTeruzzi, A., Béranger, K., Colella, S., Crise,
eling the global ocean iron cycle. Glob. BioA., 2012. Seasonal and inter-annual variabilgeochem. Cy. 18, GB1002.
ity of plankton chlorophyll and primary production in the mediterranean sea: a modelling Polimene, L., Allen, J. I., Zavatarelli, M., Jun
2006. Model of interactions between dissolved
approach. Biogeosciences 9 (1), 217–233.
organic carbon and bacteria in marine systems.
Lefevre, N., Watson, A. J., 1999. Modeling the
Aquat. Microb. Ecol. 43 (2), 127–138.
geochemical cycle of iron in the oceans and
its impact on atmospheric CO2 concentrations. Reinart, A., Arst, H., Blanco-Sequeiros, A., Herlevi, A., 1998. Relation between underwater
Global Biogeochem Cy 13, 727–736.
irradiance and quantum irradiance in depenLegendre, L., Ackley, S. F., Dieckmann, G. S.,
dence on water transparency at different depths
Gulliksen, B., Horner, R., Hoshiai, T., Melin the water bodies. J. Geophys. Res. 103 (C4),
nikov, I. A., Reeburgh, W. S., Spindler, M.,
7759–7752.
Sullivan, C. W., 1992. Ecology of sea ice biota:
2. global significance. Polar Biology 12, 429– Ruardij, P., Van Raaphorst, W., 1995. Benthic nutrient regeneration in the ERSEM ecosystem
444.
model of the North Sea. J. Sea Res. 33 (3-4),
Lengaigne, M., Menkes, C., Aumont, O.,
453–483.
Gorgues, T., Bopp, L., Madec, J.-M. A. G.,
2007. Bio-physical feedbacks on the tropical Sakshaug, E., Bricaud, A., Dandonneau, Y.,
Falkowski, P. G., Kiefer, D. A., Legendre, L.,
pacific climate in a coupled general circulation
Morel, A., Parslow, J., Takahashi, M., 1997.
model. Clim. Dyn. 28, 503–516.
Parameters of photosynthesis: definitions, theLizotte, M. P., Sullivan, C. W., 1991.
ory and interpretation of results. J. Plankt. Res.
Photosynthesis-irradiance relationships in
19, 1637–1670.
microalgae associated with Antarctic pack ice:
evidence for in situ activity. Mar. Ecol. Prog. Stoderegger, K., Herndl, G., 1998. Production
and release of bacterial capsular material and
Ser. 71, 175–184.
its subsequent utilization by marine bacterioMitra, A., Flynn, K. J., 2005. Predator-prey interplankton. Limnol. Oceanogr. 43 (5), 877–884.
actions: is ’ecological stoichiometry’ sufficient
when good food goes bad? J. Plankt. Res. 27, Sunda, W. G., Huntsman, S. A., 1997. Interrelated
influence of iron, light and cell size on marine
393–399.
phytoplankton growth. Nature 390, 389–392.
Morel, A., 1988. Optical modeling of the upper ocean in relation to its biogenous matter Tedesco, L., Vichi, M., 2010. BFM-SI: A new implementation of the biogeochemical flux model
content (case i waters). J. Geophys. Res. 93,
in sea ice. Research Papers 81, CMCC.
10,749–10,768.

103

Bibliography
Tedesco, L., Vichi, M., 2014. Sea ice biogeochemistry: A guide for modellers. PLoS ONE
9(2) (e89217).
Tedesco, L., Vichi, M., Haapala, J., Stipa, T.,
2009. An enhanced sea-ice thermodynamic
model applied to the Baltic Sea. Boreal Environmental Research 14, 68–80.

BFM Report Series 2, Bologna, Italy.
URL http://bfm-community.eu
Vichi, M., Masina, S., 11 2009. Skill assessment
of the PELAGOS global ocean biogeochemistry model over the period 1980-2000. Biogeosciences 6 (11), 2333–2353.

Vichi, M., Masina, S., Navarra, A., 2007a. A genTedesco, L., Vichi, M., Haapala, J., Stipa, T.,
eralized model of pelagic biogeochemistry for
2010. A dynamic biologically-active layer for
the global ocean ecosystem. Part II: numerical
numerical studies of the sea ice ecosystem.
simulations. J. Mar. Sys. 64, 110–134.
Ocean Modelling 35, 89–104.
Vichi, M., Pinardi, N., Masina, S., 2007b. A generalized model of pelagic biogeochemistry for
Tedesco, L., Vichi, M., Thomas, D. N., 2012. Prothe global ocean ecosystem. Part I: theory. J.
cess studies on the ecological coupling between
Mar. Sys. 64, 89–109.
sea ice algae and phytoplankton. Ecol. Mod.
226, 120–138.
Vichi, M., Ruardij, P., Baretta, J. W., 2004. Link
or sink: a modelling interpretation of the open
Thomas, D., Dieckmann, G., 2002. Antartic sea
Baltic biogeochemistry. Biogeosciences 1, 79–
ice - a habitat for extremophiles. Science 295,
100.
5555, 641–644.
Tsiaras, K. P., Kourafalou, V. H., Davidov, A., Wolf-Gladrow, D. A., Zeebe, R. E., Klaas, C.,
Staneva, J., 2008. A three-dimensional coupled
Körtzinger, A., Dickson, A. G., 2007. Total
model of the western black sea plankton dyalkalinity: The explicit conservative expresnamics: Seasonal variability and comparison
sion and its application to biogeochemical proto seawifs data. J. Geophys. Res. 113 (C7),
cesses. Marine chemistry 106 (1), 287–300.
C07007.
Worden, A. Z., Nolan, J. K., Palenik, B., 2004.
Varela, R., Cruzado, A., Gabaldon, J., 1995. ModAssessing the dynamics and ecology of maelling primary production in the North Sea
rine picophytoplankton: The importance of the
using the European Regional Seas Ecosystem
eukaryotic component. Limnol. Oceanogr. 49,
Model. J. Sea Res. 33 (3-4), 337–361.
168–179.
Vichi, M., Lovato, T., Gutierrez Mlot, E., McK- Zeebe, R. E., Wolf-Gladrow, D. A., 2001. CO2
iver, W., April 2015. Coupling BFM with ocean
in Seawater: Equilibrium, Kinetics, Isotopes.
models: the NEMO model (Nucleus for the EuVol. 65 of Oceanography Book Series. Elsevier,
ropean Modelling of the Ocean). Release 1.0.
Amsterdam.

104



---

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.4
Linearized                      : No
Page Count                      : 104
Page Mode                       : UseOutlines
XMP Toolkit                     : XMP toolkit 2.9.1-13, framework 1.6
About                           : uuid:a13e2a8b-88ad-11f0-0000-b599c202163e
Producer                        : dvips + GPL Ghostscript 9.07
Keywords                        : ()
Modify Date                     : 2015:09:01 12:06:15+02:00
Create Date                     : 2015:09:01 12:06:15+02:00
Creator Tool                    : LaTeX with hyperref package
Document ID                     : uuid:a13e2a8b-88ad-11f0-0000-b599c202163e
Format                          : application/pdf
Title                           : ()
Creator                         : ()
Description                     : ()
Subject                         : 
Author                          : 

EXIF Metadata provided by EXIF.tools