# Theoretical Analysis and Optimization of Nonlinear ODE Systems for Marine Ecosystem Models

Anna Heinle, Thomas Slawig

► To cite this version:

Anna Heinle, Thomas Slawig. Theoretical Analysis and Optimization of Nonlinear ODE Systems for Marine Ecosystem Models. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.501-510, 10.1007/978-3-642-36062-6_50. hal-01347576

HAL Id: hal-01347576

https://inria.hal.science/hal-01347576

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Theoretical analysis and optimization of
nonlinear ODE systems for marine ecosystem
models

Anna Heinle and Thomas Slawig *

Institute for Computer Science
Cluster "The Future Ocean", Christian-Albrechts Universität zu Kiel, Germany
{ahe,ts}@informatik.uni-kiel.de

**Abstract.** We present the investigation of a biogeochemical marine ecosystem model used as part of the climate change research focusing on the enhanced carbon dioxide concentration in the atmosphere. Numerical parameter optimization has been performed to improve representation of observational data using data assimilation techniques. Several local minima were found but no global optimum could be identified. To detect the actual capability of the model in simulating natural systems, a theoretical analysis of the model equations is conducted. Here, basic properties such as continuity and positivity of the model equations are investigated.

**Keywords:** Climate models, Marine ecosystem models, Parameter optimization, Ordinary differential equations

# 1 Introduction

An important part of climate change research is the investigation of biogeochemical processes occurring in the oceans. The Earth's carbon cycle, one of the main climate drivers, is highly dependent on the marine ecosystem and its interactions due to primary carbon producers such as phytoplankton. For this reason, scenarios of the climate's future are commonly created from numerical models, including a submodel to simulate the marine ecosystem. These submodels range from conceptually simple models, like the FDM model by Fasham [1], to highly complex models simulating numerous components of the marine ecosystem such as different types of plankton or multiple nutrients (see e.g. [2-4]).

In this work, we consider a model of mid-complexity used at GEOMAR, Kiel studying the CO₂ uptake of the ocean. Four components of the marine ecosystem, namely nitrogen (N), phytoplankton (P), zooplankton (Z) and detritus (D) are simulated in this model. The model, hereafter called NPZD model, depicts the main interactions and feedbacks of the marine ecosystem. However, its feasibility to reproduce real observed data is limited (see e.g. [5, 6]). Currently, the reason for this fact is unclear and is the motivation for this study.

<footnote>* In cooperation with Andreas Oschlies, GEOMAR, Kiel, Germany.</footnote>
We introduce two methods used. Numerical parameter optimization is conducted to assess the influence of different parameter combinations on the model outcome and to detect whether optimal parameters exist reproducing a given observational data set. Further, the theoretical framework of the model is analysed. The knowledge of basic properties of the model as well as the demonstration of particular dependencies of model parameters may yield a better understanding and assessment of the model dynamics.

The paper is structured as follows: In Section 2, we introduce the ecosystem model considered. Biological definitions as well as a mathematical formulation of the model are provided. In Section 3 we show an extract of the numerical experiments conducted focusing on the variability of model outcomes. Theoretical findings regarding the mathematical framework of the NPZD model are subsequently presented in Section 4. A final discussion and an outlook of future work conclude the paper in Section 5.

## 2 Model

We study the behaviour of a marine ecosystem model of NPZD type developed by Oschlies and Garçon [7]. This model simulates the concentration (in mmol Nm⁻³) of four geobiochemical components (also named *tracers* in the following) in the atlantic ocean, namely dissolved inorganic nitrogen (N), phytoplankton (P), zooplankton (Z) and detritus (D). The interactions of the four tracers among each other control the dynamics of the model (see below and Fig. 1).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Scheme of the coupling between the model variables N, P, Z and D. Arrows indicate a nitrogen flux from one to another component.</figcaption></figure>

Aside from light, phytoplankton needs nutrients to photosynthesize and to grow. In the model, these nutrients are represented by the nitrogen component N. Zooplankton graze on phytoplankton and fecal particles of zooplankton as well as other organic, sinking material are summarized in the component detritus. A part of the detritus sinks down to the bottom of the system, another part is recycled by bacteria and comes back to the nitrogen component.
## 2.1 Mathematical formulation

The model is directly defined by the interactions of the four tracers N, P, Z and D, in the following denoted by $\mathbf{y} = (y_i)_{i=1,...,4} = (N, P, Z, D)$.

A system of coupled, nonlinear ordinary differential equations (ODEs), in detail given by

$$ 
\left\{ 
\begin{aligned}
\frac{\partial y_1}{\partial t} &= -J(\mu, U)y_2 + \Phi_m^Z y_3 + \gamma_m y_4, \\
\frac{\partial y_2}{\partial t} &= (J(\mu, U) - \Phi_m^P)y_2 - G(\epsilon, g)y_3, \\
\frac{\partial y_3}{\partial t} &= (\beta G(\epsilon, g) - \Phi_m^Z - \Phi_Z^*) y_3, \\
\frac{\partial y_4}{\partial t} &= \Phi_m^P y_2 + ((1-\beta)G(\epsilon, g) + \Phi_Z^*) y_3 - (\gamma_m + w_s)y_4
\end{aligned}
\right. , \quad (1)
$$

describes these interactions. $J$ and $G$ are nonlinear functions representing the growth of phytoplankton respectively the grazing of zooplankton on phytoplankton. Parameters appearing in the equations as well as in functions $J$ and $G$ are nonnegative and signify for example growth and mortality rates. Table 1 gives a summary of all parameters. Characteristic, environmental conditions are taken into account by an additional extrinsic physical forcing.

We want to note, that due to our focus on the model specific representation of the interactions of the four tracers this forcing is not further regarded in this work.

<table><caption>**Table 1.** Model parameters. Units and definitions of the parameters are given together with their value resp. the biological, significant range. Parameters with an index in the first column are included in the optimization processes.</caption><thead><tr><th>Index</th><th>Symbol</th><th>Value</th><th>Unit</th><th>Definition</th></tr></thead><tbody><tr><td></td><td>C<sub>ref</sub></td><td>1.066</td><td>1</td><td>Growth coefficient of phytoplankton</td></tr><tr><td></td><td>c</td><td>1</td><td>°C</td><td>Growth coefficient of phytoplankton</td></tr><tr><td></td><td>k<sub>water</sub></td><td>25</td><td>m<sup>-1</sup></td><td>PAR extinction length</td></tr><tr><td></td><td>f<sub>PAR</sub></td><td>0.43</td><td>1</td><td>PAR fraction of insolation</td></tr><tr><td>1</td><td>β</td><td>[0,1]</td><td>1</td><td>Assimilation efficiency of zooplankton</td></tr><tr><td>2</td><td>ν<sub>m</sub></td><td>R<sub>+</sub></td><td>d<sup>-1</sup></td><td>Phytoplankton growth rate</td></tr><tr><td>3</td><td>α</td><td>R<sub>+</sub></td><td>m<sup>2</sup> W<sup>-1</sup>d<sup>-1</sup></td><td>Slope of photosynthesis vs light intensity</td></tr><tr><td>4</td><td>Φ<sup>Z</sup><sub>m</sub></td><td>[0,1]</td><td>d<sup>-1</sup></td><td>Zooplankton linear loss rate</td></tr><tr><td>5</td><td>k<sub>P</sub></td><td>[0, 1]</td><td>m<sup>2</sup> (mmol N)<sup>-1</sup></td><td>Light attenuation by phytoplankton</td></tr><tr><td>6</td><td>ε</td><td>R<sub>+</sub></td><td>m<sup>6</sup> (mmol N)<sup>-2</sup> d<sup>-1</sup></td><td>Grazing encounter rate</td></tr><tr><td>7</td><td>g</td><td>R<sub>+</sub></td><td>d<sup>-1</sup></td><td>Maximum grazing rate</td></tr><tr><td>8</td><td>Φ<sup>P</sup><sub>m</sub></td><td>[0,1]</td><td>d<sup>-1</sup></td><td>Phytoplankton linear mortality</td></tr><tr><td>9</td><td>Φ<sup>Z</sup><sub>Z</sub></td><td>R<sub>+</sub></td><td>m<sup>3</sup> (mmol N)<sup>-1</sup> d<sup>-1</sup></td><td>Zooplankton quadratic mortality</td></tr><tr><td>10</td><td>γ</td><td>[0,1]</td><td>d<sup>-1</sup></td><td>Detritus remineralization rate</td></tr><tr><td>11</td><td>k<sub>N</sub></td><td>R<sub>+</sub></td><td>mmol N m<sup>-3</sup></td><td>Half saturation for NO3 uptake</td></tr><tr><td>12</td><td>w<sub>s</sub></td><td>[0,1]</td><td>m d<sup>-1</sup></td><td>Detritus sinking</td></tr></tbody></table>
## Light and nutrient limited growth rate of phytoplankton

Biologically, the growth of phytoplankton, in the model represented by $J$, is limited by two factors, light and nutrients. Here, this limitation is modeled using Liebig's law of the minimum (see Eq. (2)). The temperature weighted maximum growth rate $V_p(t) = \mu_m C^{T(t)}$ is multiplied by the minimum of two functions, function $\mu$, describing the light limitation (Eq. (3)) and function $U$, mirroring the nutrient limitation (Eq. (4)).

$$J(\mu(I, y_2, t), U(y_1, t)) = V_p(t) \min\{\mu(I, y_2, t), U(y_1, t)\}. \quad (2)$$

The function $\mu$, given by

$$\mu(I, y_2, t) = \frac{1}{k_{to}(t)z} \ln \left( \frac{I_0(t) + \sqrt{V_p(t)^2 + I_0(t)^2}}{I_z(t) + \sqrt{V_p(t)^2 + I_z(t)^2}} \right), \quad (3)$$

$$I_0(t) = \alpha I_{in}(t),$$

$$I_z(t) = I_0(t)e^{-k_{to}z(t)},$$

$$k_{to}(t) = k_w(t) + k_P(t)y_2,$$

is based on the Smith's curve [8] and describes the relationship between photosynthesis and light for phytoplankton. The motivation for this function is to enable the integration of variable insolation, for example due to seasonal influences or in case a daily cycle is to be simulated.

The nutrient limitation $U$ is described by a so called Holling type II function,

$$U(y_1, t) = \frac{y_1}{k_N + y_1}. \quad (4)$$

Holling type functions describe the reproduction of a consumer as a function of food density, here as a function of nutrients $y_1$. Obviously, $U$ is monotonically increasing and ranges in the interval $[0, 1]$ for $k_N > 0$ and $y_1 \in \mathbb{R}_+^0$.

## Zooplankton growth rate

The growth of zooplankton is dependent on the phytoplankton availability and is given by a Holling type III function,

$$G(\epsilon, g, y_2) = \frac{g\epsilon y_2^2}{g + \epsilon y_2^2}. \quad (5)$$

Such as $U$, $G$ is monotonically increasing, but ranges in the interval $[0, g]$ for $g, \epsilon > 0$ and $y_2 \in \mathbb{R}_+^0$.
## 3 Numerical Optimization

We show an extract of the experiments conducted in the course of the parameter optimization of the ecosystem model introduced in Section 2. Two observational data sets are used to optimize the parameters of the model by data assimilation techniques (see below). In both cases, a cost function of least-squares type,

$$ \min_{u \in \mathbb{R}_+^k} J(y(u)) := ||y(u) - y_d||^2 \text{ s.t. } lb \le u \le ub, \quad (6) $$

is applied. Reference data is denoted by $y_d$, the model output is given by $y(u)$ and $u \in \mathbb{R}_+^k$ is the parameter vector to be optimized, componentwise within the constraints $lb$ and $ub$.

### 3.1 Data

The following two data sets are used for parameter optimization of the NPZD model.

**D1** The Bermuda Atlantic Time Series (BATS): Vertical profiles obtained during the US JGOFS project. This data is frequently used within optimization approaches of marine ecosystem models, among others the one dimensional version of the NPZD model (see e.g. [9, 6, 12]). For this study, each profile is vertically averaged from the surface to the respective “mixed layer depth” which is fixed by a temperature decline of more than 0.5 °C (for details see [10]).

**D2** Indoor-mesocosm data obtained during the AQUASHIFT project: A 30 day mesocosm experiment performed at GEOMAR, Kiel investigating the impact of ocean warming on the phytoplankton spring bloom in the North Atlantic (for details see [11]).

Since parameter optimization with respect to D1 is subject of [12], we here focus on the optimization using data set D2.

### 3.2 Results with respect to D2

In the following we show a cut-out of the numerical optimization experiments. Note, that we show individual examples instead of statistics to highlight specific issues. Table 2 and Figure 2 demonstrate the impact of varying upper bounds of the parameters (E1), Table 3 and Figure 3 provide the results of an initial value experiment (E2).

All results shown are obtained operating an optimization algorithm of SQP (sequential quadratic programming) type. Initial values are taken from literature as well as chosen randomly. The optimization runs include the indexed parameters presented in Table 1.
<table><caption>**Table 2.** Variation of the upper bounds. Initial values and upper bounds of the parameters according to [6] (left), resp. modified (right).</caption><thead><tr><th></th><th>lb</th><th>ub</th><th>u<sup>ini</sup></th><th>u<sup>*</sup></th><th>ub<sub>1</sub></th><th>ũ<sup>ini</sup></th><th>ũ<sup>*</sup></th></tr></thead><tbody><tr><td>Φ<sub>P</sub></td><td>0.0001</td><td>1</td><td>0.001</td><td>0.0021</td><td>1</td><td>0.001</td><td>0.0205</td></tr><tr><td>Φ<sup>*</sup><sub>Z</sub></td><td>0.0001</td><td>10</td><td>0.202</td><td>0.1631</td><td>1</td><td>0.202</td><td>0.4179</td></tr><tr><td>γ</td><td>0.0001</td><td>1</td><td>0.092</td><td>0.0044</td><td>1</td><td>0.092</td><td>0.9051</td></tr><tr><td>β</td><td>0.0001</td><td>1</td><td>1.000</td><td>0.9484</td><td>1</td><td>0.8</td><td>0.9365</td></tr><tr><td>Φ<sub>Z</sub></td><td>0.0001</td><td>1</td><td>0.025</td><td>0.0087</td><td>1</td><td>0.025</td><td>0.4074</td></tr><tr><td>g</td><td>0.0001</td><td>100</td><td>20.00</td><td>20.457</td><td>1</td><td>0.8</td><td>0.9761</td></tr><tr><td>ε</td><td>0.0001</td><td>100</td><td>5.446</td><td>4.6083</td><td>1</td><td>0.8</td><td>0.3319</td></tr><tr><td>ν<sub>m</sub></td><td>0.0001</td><td>10</td><td>1.076</td><td>1.1482</td><td>1</td><td>0.8</td><td>0.1696</td></tr><tr><td>k<sub>N</sub></td><td>0.0001</td><td>10</td><td>1.827</td><td>1.1982</td><td>1</td><td>0.8</td><td>0.5867</td></tr><tr><td>w<sub>s</sub></td><td>0.0001</td><td>1</td><td>0.230</td><td>0.0085</td><td>1</td><td>0.8</td><td>0.3761</td></tr><tr><td>α</td><td>0.0001</td><td>100</td><td>0.107</td><td>0.0217</td><td>1</td><td>0.107</td><td>0.2462</td></tr><tr><td>k<sub>P</sub></td><td>0.0001</td><td>1</td><td>0.026</td><td>0.8503</td><td>1</td><td>0.026</td><td>0.0130</td></tr><tr><td>Cost</td><td></td><td></td><td></td><td>3.4240</td><td></td><td></td><td>3.8091</td></tr></tbody></table>

**E1** In the first experiment, optimal parameter values found by Rückelt et al. [6] are used to initialize the model. Table 2 shows the parameters to be optimized, their initial values $u^{ini}$, optimized values $u^*$ and the bounds $lb$ and $ub$ for two examples. Replacing the upper bounds as given in [6] by a vector of ones$^1$, the cost function can be decreased less ($J = 3.8091$) in contrast to final costs using the original setting ($J = 3.4240$). One might assume that narrowing the constraints may yield to worse optimization results. However, the values do not differ so much, and, as visible in Figure 2, both parameter sets result in an unsatisfying model outcome.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Simulations of N, P, Z and D for a 30 day period, according to the length of the mesocosm study, using parameter values $u^*$ (left) and $\tilde{u}^*$ (right) as presented in Table 2. Circles depict data points.</figcaption></figure>

**E2** The second experiment addresses the influence of initial values on the optimized parameter vector. Table 3 displays the results of four, randomly initialized optimization runs, where $u_{1,2}^{ini} \in [0, 1]^{12}$ and $u_{3,4}^{ini}$ are within the constraints (see Tab. 3). While the optimization started at $u_1^{ini}$ reaches costs $J$ near 1, the optimization process initialized with $u_2^{ini}$ gets stuck, just reaching a cost function of $J \approx 6$. Similar results are obtained for $u_3^{ini}$ and $u_4^{ini}$ indicating the independence of the optimization on the initial values. Further,

<footnote n="1">¹ Consequently, the initial values of some parameters are modified as well.</footnote>
although the cost function reaches values in the same ranges, the outcomes as well as the final parameter sets $u_i^*$, $i = 1, \dots, 4$ differ noticeably and we deduce the existence of multiple distinctive, local minima (see Tab. 3 and Fig. 3).

<table><caption>**Table 3.** Impact of different, random initial values $u_i^{ini}$. Upper and lower bounds are fix ($lb = (0.0001)^{12}$ and $ub_i = 1.0$ for $i = 1, 3, 4, 5, 12$, $ub_j = 10.0$, $j \neq i$). $u^*$ are the optimized values.</caption><thead><tr><th></th><th>u<sub>1</sub><sup>ini</sup></th><th>u<sub>1</sub><sup>*</sup></th><th>u<sub>2</sub><sup>ini</sup></th><th>u<sub>2</sub><sup>*</sup></th><th>u<sub>3</sub><sup>ini</sup></th><th>u<sub>3</sub><sup>*</sup></th><th>u<sub>4</sub><sup>ini</sup></th><th>u<sub>4</sub><sup>*</sup></th></tr></thead><tbody><tr><td>Φ<sub>P</sub></td><td>0.8147</td><td>0.0585</td><td>0.6787</td><td>0.5547</td><td>0.1146</td><td>0.0036</td><td>0.4551</td><td>0.0005</td></tr><tr><td>Φ<sub>Z</sub><sup>*</sup></td><td>0.9058</td><td>0.0660</td><td>0.7577</td><td>1.7077</td><td>3.5730</td><td>8.8814</td><td>3.7551</td><td>3.3658</td></tr><tr><td>γ</td><td>0.1270</td><td>0.1883</td><td>0.7431</td><td>0.9756</td><td>0.0947</td><td>0.6833</td><td>0.0095</td><td>0.1377</td></tr><tr><td>β</td><td>0.9134</td><td>0.7945</td><td>0.3922</td><td>0.4261</td><td>0.0942</td><td>0.2343</td><td>0.2697</td><td>0.1882</td></tr><tr><td>Φ<sub>Z</sub></td><td>0.6324</td><td>0.0035</td><td>0.6555</td><td>0.2604</td><td>0.0476</td><td>0.5000</td><td>0.1297</td><td>0.9283</td></tr><tr><td>g</td><td>0.0975</td><td>0.0361</td><td>0.1712</td><td>3.4484</td><td>7.0978</td><td>3.4276</td><td>6.3543</td><td>7.9393</td></tr><tr><td>ε</td><td>0.2785</td><td>0.7821</td><td>0.7060</td><td>1.1090</td><td>0.9482</td><td>1.3885</td><td>2.4897</td><td>0.2272</td></tr><tr><td>ν<sub>m</sub></td><td>0.5469</td><td>0.1146</td><td>0.0318</td><td>0.7970</td><td>0.6744</td><td>0.0297</td><td>4.2283</td><td>4.6502</td></tr><tr><td>k<sub>N</sub></td><td>0.9575</td><td>3.0715</td><td>0.2769</td><td>1.8394</td><td>6.6350</td><td>6.5590</td><td>1.3252</td><td>2.2505</td></tr><tr><td>w<sub>s</sub></td><td>0.9649</td><td>0.0458</td><td>0.0462</td><td>0.3040</td><td>9.7396</td><td>9.2367</td><td>4.8159</td><td>5.7531</td></tr><tr><td>α</td><td>0.1576</td><td>9.9534</td><td>0.0971</td><td>0.1138</td><td>5.8376</td><td>5.8376</td><td>2.1038</td><td>2.2880</td></tr><tr><td>k<sub>P</sub></td><td>0.9706</td><td>0.6453</td><td>0.8235</td><td>0.9280</td><td>0.0162</td><td>0.0162</td><td>0.5233</td><td>0.5997</td></tr><tr><td>Cost</td><td></td><td>1.109</td><td></td><td>5.925</td><td></td><td>1.390</td><td></td><td>5.670</td></tr></tbody></table>

These two experiments demonstrate the general challenge of finding an optimal parameter set in the context of ecosystem modeling. The experiments performed for the NPZD model indicate the existence of numerous local minima which make the numerical detection of a global minimum - as far as it exists - virtually impossible. Thus, we proceed with a second approach to gain information on the feasibility of the model.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Simulations created by $u_1^*$ (top left), $u_2^*$ (top right), $u_3^*$ (bottom left) and $u_4^*$ (bottom right). Notations are the same as in Figure 1.</figcaption></figure>
# 4 Theoretical analysis

In the following, we consider the initial value problem (IVP)

$$ \mathbf{y}' = \mathbf{f}(t, \mathbf{y}), \mathbf{y}(t_0) = \mathbf{y}_0 > 0 \text{ and } t \in I = [t_0, t_0 + a], \quad (7) $$

where $\mathbf{f}(t, \mathbf{y})$ is given by the right-hand side of

$$ \left\{ \begin{aligned} \frac{\partial y_1}{\partial t} &= -J(U)y_2 + \Phi_m^Z y_3 + \gamma_m y_4, \\ \frac{\partial y_2}{\partial t} &= (J(U) - \Phi_m^P)y_2 - G(\epsilon, g)y_3, \\ \frac{\partial y_3}{\partial t} &= (\beta G(\epsilon, g) - \Phi_m^Z - \Phi_Z^* y_3)y_3, \\ \frac{\partial y_4}{\partial t} &= \Phi_m^P y_2 + ((1 - \beta)G(\epsilon, g) + \Phi_Z^* y_3)y_3 - \gamma_m y_4 \end{aligned} \right. . \quad (8) $$

In contrast to system (1), we here assume that the ODE system is mass conserving (equivalent to $w_s = 0$ in system (1)) and the limitation of phytoplankton growth due to light is ignored.

To reveal fundamental characteristics of the model dynamics, basic properties of the model equations have to be checked at first. From now on, the differential equations $\frac{\partial y_i}{\partial t}$ are denoted by $f_i(\mathbf{y})$, $i = 1, \dots, 4$ and the model parameters are summarized in the vector $\mathbf{u} = (u_i)_{i=1,\dots,12}$.

**Remark 1** The equations $f_i$, $i = 1, \dots, 4$ in (8) are continuous in $\mathbf{y}$ and $\mathbf{u}$ for positive $\mathbf{y}$ resp. $\mathbf{u}$.

This is obvious since the equations are compositions of continuous functions in both, $\mathbf{y}$ and $\mathbf{u}$, that especially holds for the functions $J$ and $G$.

Next, we consider the range of $\mathbf{f}$. Concentrations, as simulated in our model, are assumed to be nonnegative. To ensure this, we check an important preliminary to investigate this property. The following definition is used.

**Definition 1** A function $f: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n$ is called *quasipositive* if for $\mathbf{y} \in \mathbb{R}^n$ with $y_i \ge 0$, $i = 1, \dots, n$, for all $k \in 1, \dots, n$, $t \ge t_0$ and $y_k = 0$ hold $f_k(t, \mathbf{y}) \ge 0$.

**Proposition 1** The equations $f_i$, $i = 1, \dots, 4$ in (8) are quasipositive given any parameter vector $\mathbf{u} \ge 0$.

*Proof.* Let $\mathbf{u} \ge 0$ and $\mathbf{y} = (y_i)_{i=1,\dots,4} \ge 0$. Then

$$ \begin{aligned} f_1(\mathbf{y}) &= \Phi_m^Z y_3 + \gamma_m y_4 \ge 0 \text{ for } y_1 = 0, \\ f_2(\mathbf{y}) &= 0 \text{ for } y_2 = 0, \\ f_3(\mathbf{y}) &= 0 \text{ for } y_3 = 0, \\ f_4(\mathbf{y}) &= \Phi_m^P y_2 + ((1 - \beta)G(\epsilon, g) + \Phi_Z^* y_3)y_3 \ge 0 \text{ for } y_4 = 0. \end{aligned} $$

□

To investigate the existence and uniqueness of solutions of an IVP, Lipschitz continuity plays a main role. We want to recall the following, well-known Lemma.
**Lemma 1** Let $\mathcal{D}$ be a domain in $\mathbb{R}^n$ and $f: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n$. If $f$ and its partial derivatives $\frac{\partial f_i}{\partial y_j}$ are continuous in $\mathcal{D}$, then $f$ is locally Lipschitz continuous in $\mathcal{D}$.

Now, for parameters in the range $[0, 1]$, the local Lipschitz continuity of functions $f_i$, $i = 1, \dots, 4$ can be shown providing a good base for further analyses regarding the uniqueness of solutions.

**Proposition 2** Let $\mathcal{D} \subset \mathbb{R}_+^4$ be a domain and $\mathbf{u} \in [0, 1]^{12}$. Then, the equations $f_i, i = 1, \dots, 4$ in (8) are locally Lipschitz continuous in $\mathbf{y}$.

*Proof.* As noted before, the functions $f_i$ are continuous for all $i$. We do not show the calculations in detail, but $\frac{\partial f_i}{\partial y_j} \in \mathcal{C}^1(\mathcal{D})$ for $i, j = 1, \dots, 4$. Hence, $f \in \mathcal{C}^1(\mathcal{D})$ and the preliminaries of corollary 1 are fulfilled. The local Lipschitz continuity of $f$ follows. $\square$

**Remark 2** The global Lipschitz continuity of $f_i, i = 1, \dots, 4$ is not necessarily given. For demonstration, we here give a short counterexample.

**Example:** Let $L > 0$ and a parameter vector $\mathbf{u} \in [0, 1]$. We distinguish 2 cases.

1. $\beta g \le \beta g^2 + 2\Phi_Z^* + \Phi_Z$.

Choose $\mathbf{y} := (y_1, \sqrt{\frac{1-g}{\epsilon}}, 1, y_4)$ and $\bar{\mathbf{y}} := (y_1, \sqrt{\frac{1-g}{\epsilon}}, 1 + \frac{L+\delta}{\Phi_Z^*}, y_4)$ for a $\delta > 0$. Then

$$ \begin{aligned} ||f_3(\mathbf{y}) - f_3(\bar{\mathbf{y}})|| &= (L + \delta)||\frac{\Phi_Z}{\Phi_Z^*} + 2 + \frac{L + \delta}{\Phi_Z^*} - \frac{\beta g(1 - g)}{\Phi_Z^*}|| \\ &> L\frac{L + \delta}{\Phi_Z^*} = L||\mathbf{y} - \bar{\mathbf{y}}||. \end{aligned} $$

2. $\beta g > \beta g^2 + 2\Phi_Z^* + \Phi_Z$.

Choose $\mathbf{y} := (y_1, \sqrt{\frac{1-g}{\epsilon}}, 1, y_4)$ and $\bar{\mathbf{y}} := (y_1, \sqrt{\frac{1-g}{\epsilon}}, 1 - \frac{L}{\Phi_Z^*}, y_4)$. Then

$$ \begin{aligned} ||f_3(\mathbf{y}) - f_3(\bar{\mathbf{y}})|| &= L||\frac{\Phi_Z}{\Phi_Z^*} + 2 + \frac{L}{\Phi_Z^*} - \frac{\beta g(1 - g)}{\Phi_Z^*}|| \\ &> L\frac{L}{\Phi_Z^*} = L||\mathbf{y} - \bar{\mathbf{y}}||. \end{aligned} $$

# 5 Conclusions & Outlook

An extract of the numerical and theoretical investigation of a marine ecosystem model is presented. Parameter optimization and numerical experiments are
conducted to get an insight of the model variability. Two examples are shown
addressing the influence of different parameter settings. We demonstrate that
the model is highly variable and numerous local minima exist, making the de-
termination of a global optimum by numerical methods virtually impossible.
The theoretical part of this study focus on basic properties of the model equa-
tions. Continuity and quasipositivity of the equations are established and con-
ditions for local Lipschitz continuity are presented. These first findings facilitate
the ongoing investigation of theoretical characteristics of the model making the
model and its application in climate research more meaningful.

References

1. Fasham M.J.R., Ducklow, H.W., McKelvie, S.M.: A nitrogen-based model of plankton dynamics in the oceanic mixed layer. J. Mar. Res. 99, 591–639 (1990)

2. Bissett, W.P., Walsh, J.J., Dieterle, D.A., Carder, K. L.: Carbon cycling in the upper waters of the Sargasso Sea: I. Numerical simulation of differential carbon and nitrogen fluxes. Deep-Sea Res. I 46, 205–269 (1999)

3. Moore, J.K., Doney, S.C., Kleypas, J.A., Glover, D.M., Fung, I.Y.: An intermediate complexity marine ecosystem model for the global domain. Deep-Sea Res. II 49, 403462. (2002)

4. Lancelot, C., Spitz, Y.H., Gypens, N., Ruddick, K., Becquevort, S., Rousseau, V., Billen, G.: Modelling diatom-Phaeocystis blooms and nutrient cycles in the Southern Bight of the North Sea: the MIRO model. Mar. Ecol. Prog. Ser. 289, 63–78 (2005)

5. Schartau, M., Oschlies, A.: Simultaneous data-based optimization of a 1d-ecosystem model at three locations in the north atlantic: Part I - method and parameter estimates. J. Mar. Res. 61, 765–793 (2003)

6. Rückelt, J., Sauerland, V., Slawig, T., Srivastav, A., Ward, B., Patvardhan, C.: Parameter optimization and uncertainty analysis in a model of oceanic CO2-uptake using a hybrid algorithm and algorithmic differentiation. Nonlinear Anal. Real, Online (2010)

7. Oschlies, A., Garçon, V.: An eddy-permitting coupled physical-biological model of the north atlantic. 1. Sensitivity to advection numerics and mixed layer physics. Global Biogeochem. Cy. 13, 135–160 (1999)

8. Jassby, A.D., Platt, T.: Mathematical formulation of the relationship between photo synthesis and light for phytoplankton. Limnol. Oceanogr. 21, 540–547 (1976)

9. Spitz, Y.H., Moisan, J.R., Abbott, M.R.: Configuring an ecosystem model using data from the Bermuda Atlantic time series (BATS). Deep-Sea Res. II 48, 1733–1768 (2001)

10. Fasham, M.J.R., Evans, G.T.: The use of optimisation techniques to model marine ecosystem dynamics at the JGOFS station at 473N 203W. Philos. T. Roy. Soc. B 348, 206–209 (1995)

11. Sommer, U., Lengfellner, K.: Climate change and the timing, magnitude, and com-
position of the phytoplankton spring bloom. Glob. Change Biol. 14, 1199–1208
(2008)

12. El Jarbi, M., Slawig, T.: Introducing Periodic Parameters in a Marine Ecosystem Model using Linear Quadratic Optimal Control, In: Proceedings of the 25th IFIP TC7 Conference on System Modeling and Optimization (2012)