# Surrogate-based optimization using an artificial neural network for a
parameter identification in a 3D marine ecosystem model

Markus Pfeil  Kiel Marine Science (KMS) - Centre for Interdisciplinary Marine
Science, Dep. of Computer Science, Kiel University, 24098 Kiel, Germany ({mpf,
ts}@informatik.uni-kiel.de).  Thomas Slawig  Kiel Marine Science (KMS) -
Centre for Interdisciplinary Marine Science, Dep. of Computer Science, Kiel
University, 24098 Kiel, Germany ({mpf, ts}@informatik.uni-kiel.de).

###### Abstract

Parameter identification for marine ecosystem models is important for the
assessment and validation of marine ecosystem models against observational
data. The surrogate-based optimization (SBO) is a computationally efficient
method to optimize complex models. SBO replaces the computationally expensive
(high-fidelity) model by a surrogate constructed from a less accurate but
computationally cheaper (low-fidelity) model in combination with an
appropriate correction approach, which improves the accuracy of the low-
fidelity model. To construct a computationally cheap low-fidelity model, we
tested three different approaches to compute an approximation of the annual
periodic solution (i.e., a steady annual cycle) of a marine ecosystem model:
firstly, a reduced number of spin-up iterations (several decades instead of
millennia), secondly, an artificial neural network (ANN) approximating the
steady annual cycle and, finally, a combination of both approaches. Except for
the low-fidelity model using only the ANN, the SBO yielded a solution close to
the target and reduced the computational effort significantly. If an ANN
approximating appropriately a marine ecosystem model is available, the SBO
using this ANN as low-fidelity model presents a promising and computational
efficient method for the validation.

##  1 Introduction

In the field of climate research, marine ecosystem models are part and parcel
of the analysis of changes in the marine ecosystem influenced by diverse
biogeochemical processes. As part of the global carbon cycle, the ocean takes
up, for example, $\textrm{CO}_{2}$ from the atmosphere and stores it. A marine
ecosystem model consists of a global circulation model coupled with a
biogeochemical model considering the interactions of the physical and
biogeochemical processes [[, cf.]]Fasham03, SarGru06, FenNeu04. In particular,
the equations and variables describing the physical processes are well known.
Conversely, there is generally no set of equations and variables to describe
the biogeochemical processes which is why many biogeochemical models exist
that differ in their complexity by the number of state variables and
parametrizations [[, see e.g.,]]KrKhOs10. Therefore, validation and assessment
of the biogeochemical models, which include a parameter optimization and a
discussion of simulation results, are necessary by an evaluation of the model
outputs against observational data [13].

Parameter identification to validate biogeochemical models, especially for
three-dimensional models, is a challenging task with a large computationally
effort. The parameter identification determines or adjusts the model
parameters to fit the model to existing measurement data or a corresponding
model output [4]. For marine ecosystem models, this corresponds to the solving
of a nonlinear optimization problem that typically requires many evaluations
of the model-data misfit function, also called cost or objective function. As
a consequence, such optimizations using conventional optimization algorithms
are very time-consuming or infeasible even on high-performance computers
because a single model evaluation already is computationally expensive [37].
To date, the search for a feasible optimization for global marine ecosystem
models has included, for instance, approaches of a statistical emulator
technique as well as a gradient-free method [25, 26, 33]. Due to the high
computational effort, the acceleration of the optimization process itself or
of the underlying simulation is still of great interest.

_Surrogate-based optimization_ (SBO) is a computationally efficient method to
optimize complex models (so-called _high-fidelity models_), whose simple model
evaluation already requires enormous computational effort [46, 15, 3, 29]. The
surrogate, which is computationally cheap and a reasonable approximation of
the high-fidelity model, replaces the high-fidelity model in the optimization
process. Consequently, the SBO estimates the optimal parameters of the high-
fidelity model by optimizing the surrogate. There are several variants to
construct the surrogate, such as function-approximation surrogates [46, 51,
52] or physics-based ones [53] using a _low-fidelity model_ that is a less
accurate approximation of the high-fidelity model. The latter approach is used
in this paper because coarser discretizations (in time and/or space) are
common in climate research [34]. The application of SBO for the parameter
identification of marine ecosystem models has already reduced computational
costs compared to conventional optimizations because the number of high-
fidelity model evaluations were much smaller [45, 43, 44].

Several strategies enable the reduction of the computational effort simulating
marine ecosystem models. The fully coupled simulation of a marine ecosystem
model to compute a steady annual cycle is computationally expensive because it
requires a long-time integration over several millennia, but a single
evaluation is already computational expensive [[, cf.]]BeDiWu08, Bryan84,
DaMcLa96, WunHei08, SibWun11, Osc06. As opposed to the fully coupled
simulation (the so-called _online_ simulation), the _offline_ simulation
neglects the impact of the biogeochemical model on the ocean circulation to
reduce the computational effort using pre-computed data of the ocean currents.
Moreover, the _transport matrix method_ (TMM) reduces the computation of the
global ocean circulation to a matrix-vector multiplication and, thus,
decouples the evaluation of the biogeochemical model from the ocean
circulation [19, 18]. This approach also lowers the computational effort with
a tolerable loss of accuracy. These two strategies are only two examples to
reduce the computational effort.

Deep Learning enables predictions of steady annual cycles for a marine
ecosystem model using the features learned from other steady annual cycles. In
recent years, Deep Learning [16, 27] based on _artificial neural networks_
(ANNs) allowed breakthroughs in many areas, such as classification, speech
recognition, computer vision or bioinformatics. The power of Deep Learning
lies in recognizing patterns in data provided to train the neural network and
predicting corresponding data based on the learned patterns. We applied Deep
Learning to predict steady annual cycles from model parameters so that the
long-time integration was not longer necessary [40].

The combination of a low-fidelity model and correction technique affects both
the accuracy and the computational effort of the surrogate that was used for
the SBO [[, cf.]]Son03, PPKOS13. The low-fidelity model is always a compromise
between the accuracy and the computational effort of the high-fidelity model
approximation. An ANN offers new opportunities to approximate a steady annual
cycle for the use as low-fidelity model. While the accuracy using only the ANN
as low-fidelity model was not sufficient for the SBO, the SBO which used a
physics-based low-fidelity model improved by the ANN showed good results. In
the present paper, we applied synthetic data generated by the high-fidelity
model to assess the feasibility of the SBO with the different low-fidelity
models. Instead of real measurement data, the use of synthetic data avoided
uncertainties, such as the required structural complexity of the model to
reconstruct real measurement data, errors in the measurements and the
performance of the optimization itself.

This paper is structured as follows: after an introduction into marine
ecosystem models including the computation of steady annual cycles using the
TMM in Section 2, Section 3 contains the description of the methods used to
train an ANN together with the resulting ANN. In Section 4, we describe the
SBO. Numerical results of the SBO using different low-fidelity models are
presented in Section 5. The paper closes with a summary and conclusions.

##  2 Model description

A marine ecosystem model describes the interplay between the ocean circulation
and the interactions and biogeochemistry among ocean biota. The marine
ecosystem is represented by a given number of ecosystem species (or tracers)
which are substances in the ocean water and subject to chemical or biochemical
reactions. In the marine ecosystem model, both the ocean circulation affects
the tracer concentrations and, vice versa, the tracers affect the ocean
circulation. A fully coupled model (_online_ model) includes both of these
interactions and, hence, simulating such a model is computationally expensive
and often limited to single model evaluations [37]. In contrast, an _offline_
model simplifies the complexity by neglecting the influence of the tracers on
the ocean circulation. As a result of this one-way coupling, we can use a pre-
computed ocean circulation for the simulation.

###  2.1 Model equations for marine ecosystems

A system of partial differential equations models the marine ecosystem. The
complexity of the marine ecosystem models varies, inter alia, by the number of
tracers, which defines the size of the system of differential equations. The
differential equations in the system are of the Lotka-Volterra or predator-
prey type [30, 56]. In the rest of this paper, we consider marine ecosystem
models using an offline model and $n_{y}\in\mathbb{N}$ tracers on a spatial
domain $\Omega\subset\mathbb{R}^{3}$ (i.e., the ocean) and a time interval
$[0,1]$ (i.e., one model year). The mapping
$y_{i}:\Omega\times[0,1]\rightarrow\mathbb{R}$,
$i\in\left\\{1,\ldots,n_{y}\right\\}$, describes the tracer concentration of
the single tracer $y_{i}$ and $\mathbf{y}:=\left(y_{i}\right)_{i=1}^{n_{y}}$
summarizes the tracer concentrations of all tracers. The system of parabolic
partial differential equations

$\displaystyle\frac{\partial y_{i}}{\partial
t}(x,t)+\left(D(x,t)+A(x,t)\right)y_{i}(x,t)$
$\displaystyle=q_{i}\left(x,t,\mathbf{y},\mathbf{u}\right),$ $\displaystyle
x\in\Omega,t$ $\displaystyle\in[0,1],$ (1) $\displaystyle\frac{\partial
y_{i}}{\partial n}(x,t)$ $\displaystyle=0,$ $\displaystyle
x\in\partial\Omega,t$ $\displaystyle\in[0,1],$ (2)

for $i=1,\ldots,n_{y}$, defines the tracer transport of a marine ecosystem
model. The linear operators $A:\Omega\times[0,1]\rightarrow\mathbb{R}$ and
$D:\Omega\times[0,1]\rightarrow\mathbb{R}$ correspond to the advection and
diffusion coming from the ocean circulation while the term
$q_{i}:\Omega\times[0,1]\rightarrow\mathbb{R}$ summarizes all biogeochemical
terms for the tracer $y_{i}$. The homogeneous Neumann boundary condition (2)
includes the normal derivative describing an absence of fluxes on the
boundary.

The spatial tracer transport in marine water depends on the ocean currents in
form of spatially discretized advection and diffusion. The advection including
a given velocity field $v:\Omega\times[0,1]\rightarrow\mathbb{R}^{3}$ is
modeled as

$\displaystyle A(x,t)y_{i}(x,t)$
$\displaystyle:=\textrm{div}\left(v(x,t)y_{i}(x,t)\right),$ $\displaystyle
x\in\Omega,t$ $\displaystyle\in[0,1]$ (3)

for $i\in\\{1,\ldots,n_{y}\\}$. The diffusion models the turbulent effects of
the ocean circulation. The molecular diffusion of the tracers themselves,
conversely, is neglected because this diffusion is much smaller than the
diffusion induced by turbulence. In ocean circulation modeling, the quite
different scales in horizontal and vertical direction requires a splitting
$D=D_{h}+D_{v}$ and an implicit treatment of the vertical part in the time
integration. Both directions of the diffusion are modeled in the second-order
form as

$\displaystyle D_{h}(x,t)y_{i}(x,t)$
$\displaystyle:=-{\textrm{div}_{h}}\left(\kappa_{h}(x,t)\nabla_{h}y_{i}(x,t)\right)$
$\displaystyle x\in\Omega,t$ $\displaystyle\in[0,1],$ (4) $\displaystyle
D_{v}(x,t)y_{i}(x,t)$ $\displaystyle:=-\frac{\partial}{\partial
z}\left(\kappa_{v}(x,t)\frac{\partial y_{i}}{\partial z}(x,t)\right),$
$\displaystyle x\in\Omega,t$ $\displaystyle\in[0,1]$ (5)

for $i\in\\{1,\ldots,n_{y}\\}$, where $\textrm{div}_{h}$ and $\nabla_{h}$
denote the horizontal divergence and gradient,
$\kappa_{h},\kappa_{v}:\Omega\times[0,1]\rightarrow\mathbb{R}$ the diffusion
coefficient fields and $z$ is the vertical coordinate. Due to the neglect of
the molecular diffusion, the diffusion coefficients are identical for all
tracers.

The biogeochemical model summarizes the biogeochemical processes modeled in
the marine ecosystem. The nonlinear function
$q_{i}:\Omega\times[0,1]\rightarrow\leavevmode\nobreak\
\mathbb{R},\left(x,t\right)\mapsto
q_{i}\left(x,t,\mathbf{y},\mathbf{u}\right)$ describes the biogeochemical
processes for the tracer $y_{i}$, $i\in\\{1,\ldots,n_{y}\\}$. More
specifically, the biogeochemical processes depend, firstly, on the variability
of the solar radiation in space and time, secondly, on the coupling to other
tracers and, thirdly, on $n_{u}\in\mathbb{N}$ model parameters
$\mathbf{u}\in\mathbb{R}^{n_{u}}$ (such as growth, loss and mortality rates or
sinking speed). The biogeochemical model
$\mathbf{q}=\left(q_{i}\right)_{i=1}^{n_{y}}$ combines the biogeochemical
processes for all tracers. In contrast to the biogeochemical model, the marine
ecosystem model, moreover, includes the effects of the ocean circulation and,
thus, comprises the whole system (1) to (6).

For the marine ecosystem model, an annual periodic solution of (1) and (2)
(i.e., a steady annual cycle) fulfills

$\displaystyle y_{i}(x,0)$ $\displaystyle=y_{i}(x,1),$ $\displaystyle x$
$\displaystyle\in\Omega,$ (6)

for $i\in\\{1,\ldots,n_{y}\\}$. For this purpose, we assume that the operators
$D,A$ and the functions $q_{i}$ also are annually periodic in time.

###  2.2 Biogeochemical models

We applied two global biogeochemical models with increased complexity, the N
and N-DOP model [23, 42]. In the following, we briefly introduce these models.
For a detailed description of the modeled processes and model equations, we
refer to [23] and [42].

The light intensity affects the biogeochemical processes in the ocean. The
light limitation function $I:\Omega\times[0,1]\rightarrow\mathbb{R}_{\geq 0}$
models the light intensity depending on the insolation. This is based on the
astronomical formula of [38] and considers the ice cover as well as the
exponential attenuation of water. According to the light intensity, the ocean
is divided into a euphotic (sun lit) zone of about $100\text{\,}\mathrm{m}$
and an aphotic zone below. A fast and dynamic turnover of phosphorus,
especially, is part of the euphotic zone, for instance via photosynthesis,
grazing or mortality. Furthermore, a part of the biological production sinks
from the euphotic zone as a particulate matter to depth and remineralizes
there according to the empirical power-law relationship [32].

Table 1: Model parameters of the biogeochemical models. Parameter | Description | Unit
---|---|---
$k_{w}$ | Attenuation coefficient of water | ${\mathrm{m}}^{-1}$
$\mu_{P}$ | Maximum growth rate | ${\mathrm{d}}^{-1}$
$K_{N}$ | Half saturation constant for $\mathrm{P}\mathrm{O}_{\mathrm{4}}$ uptake | $\mathrm{mmol}\text{\,}\mathrm{P}\text{\,}{\mathrm{m}}^{-3}$
$K_{I}$ | Light intensity compensation | $\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$
$\sigma_{\text{DOP}}$ | Fraction of phytoplankton losses assigned to DOP | $\mathrm{1}$
$\lambda^{\prime}_{\text{DOP}}$ | Decay rate | ${\mathrm{yr}}^{-1}$
$b$ | Implicit representation of sinking speed | $\mathrm{1}$

The N model contains only one tracer describing phosphate ($\textrm{PO}_{4}$)
as inorganic nutrients [[, cf.]]BacMai90, KrKhOs10. We denote this model as N
model (N for nutrients), i.e., $\mathbf{y}=(\mathbf{y}_{\textrm{N}})$. The
$n_{u}=5$ model parameters
$\mathbf{u}=\left(k_{w},\mu_{P},K_{N},K_{I},b\right)$ (see Table 1) control
the biogeochemical processes. The phytoplankton production (or biological
uptake) is defined as

$\displaystyle f_{P}:\Omega\times[0,1]\rightarrow\mathbb{R},f_{P}(x,t)$
$\displaystyle=\mu_{P}y_{P}^{*}\frac{I(x,t)}{K_{I}+I(x,t)}\frac{\mathbf{y}_{N}(x,t)}{K_{N}+\mathbf{y}_{N}(x,t)}$
(7)

and depends on available nutrients and light. Moreover, the uptake of
nutrients by phytoplankton is limited using a half saturation function, and
applies a model parameter $\mu_{P}$ for the maximum production rate as well as
a prescribed concentration of phytoplankton $y_{P}^{*}=0.0028$
$\mathrm{mmol}\text{\,}\mathrm{P}\text{\,}{\mathrm{m}}^{-3}$.

Additional to nutrients (N), the N-DOP model includes dissolved organic
phosphorus (DOP) [[, cf.]]BacMai91, PaFoBo05, KrKhOs10, i.e.,
$\mathbf{y}=(\mathbf{y}_{\textrm{N}},\mathbf{y}_{\textrm{DOP}})$. Using the
same phytoplankton production (7) as the N model, the N-DOP contains $n_{u}=7$
model parameters
$\mathbf{u}=\left(k_{w},\mu_{P},K_{N},K_{I},\sigma_{\text{DOP}},\lambda^{\prime}_{\text{DOP}},b\right)$
(see Table 1) to describe the internal processes.

###  2.3 Transport matrix method

The _transport matrix method_ (TMM) efficiently approximates the tracer
transport by matrix-vector multiplications [19, 18]. Instead of implementing
directly a discretization scheme for the advection and diffusion operators $A$
and $D$, the TMM approximates the ocean circulation by matrices because the
application of the two operators on a spatially discretized tracer vector is
linear and, therefore, the discretized advection-diffusion equation can be
written as a linear matrix equation. For this purpose, the TMM computes and
stores the matrices obtained by applying the discretized operators on a
discrete tracer vector. Consequently, the matrices contain the transport of
all parameterized processes represented in the underlying ocean circulation
model [19].

The TMM reduces a time step of the simulation of a marine ecosystem model to
matrix-vector multiplications and an evaluation of the biogeochemical model.
We assume that, firstly, the grid with $n_{x}\in\mathbb{N}$ grid points
$\left(x_{k}\right)_{k=1}^{n_{x}}$ is a spatial discretization of the domain
$\Omega$ (i.e., the ocean) and, secondly, the time steps
$t_{0},\ldots,t_{n_{t}}\in[0,1]$, $n_{t}\in\mathbb{N}$, specified by

$\displaystyle t_{j}$ $\displaystyle:=j\Delta t,$ $\displaystyle j$
$\displaystyle=0,\ldots,n_{t},$ $\displaystyle\Delta t$
$\displaystyle:=\frac{1}{n_{t}},$

define an equidistant grid of the time interval $[0,1]$ (i.e., one model
year). For the time instant $t_{j}$, $j\in\\{0,\ldots,n_{t}-1\\}$, the vector
$\mathbf{y}_{j}:=\left(\mathbf{y}_{ji}\right)_{i=1}^{n_{y}}\in\mathbb{R}^{n_{y}n_{x}}$
combines the numerical approximations for all tracers at time instant $t_{j}$
using a reasonable concatenation whereby vector

$\displaystyle\mathbf{y}_{ji}$
$\displaystyle\approx\left(y_{i}\left(t_{j},x_{k}\right)\right)_{k=1}^{n_{x}}\in\mathbb{R}^{n_{x}}$

denotes the approximation of the spatially discrete tracer $y_{i}$,
$i\in\\{1,\ldots,n_{y}\\}$. Analogously,

$\displaystyle\mathbf{q}_{ji}$
$\displaystyle\approx\left(q_{i}\left(x_{k},t_{j},\mathbf{y_{j}},\mathbf{u}\right)\right)_{k=1}^{n_{x}}\in\mathbb{R}^{n_{x}}$

contains the spatially discretized biogeochemical term of tracer $y_{i}$,
$i\in\\{1,\ldots,n_{y}\\}$, at time instant $t_{j}$,
$j\in\\{0,\ldots,n_{t}-1\\}$, and
$\mathbf{q}_{j}:=\left(\mathbf{q}_{ji}\right)_{i=1}^{n_{y}}$ summarizes this
for all tracers at time instant $t_{j}$. In equation (1), we discretize the
advection and the horizontal diffusion explicitly and the vertical diffusion
implicitly. Using the spatially discrete counterparts
$\mathbf{A}_{j},\mathbf{D}_{j}^{h}$ and $\mathbf{D}_{j}^{v}$ of the operators
$A,D_{h}$ and $D_{v}$ at time instant $t_{j}$, $j\in\\{0,\ldots,n_{t}-1\\}$,
the semi-implicit Euler scheme of (1) results in a time-stepping

$\displaystyle\mathbf{y}_{j+1}$ $\displaystyle=\left(\mathbf{I}+\Delta
t\mathbf{A}_{j}+\Delta t\mathbf{D}_{j}^{h}\right)\mathbf{y}_{j}+\Delta
t\mathbf{D}_{j}^{v}\mathbf{y}_{j+1}+\Delta
t\mathbf{q}_{j}\left(\mathbf{y}_{j},\mathbf{u}\right),$ $\displaystyle j$
$\displaystyle=0,\ldots,n_{t}-1$

with the identity matrix $\mathbf{I}\in\mathbb{R}^{n_{x}\times n_{x}}$.
Defining the explicit and implicit transport matrices

$\displaystyle\mathbf{T}_{j}^{\text{exp}}$ $\displaystyle:=\mathbf{I}+\Delta
t\mathbf{A}_{j}+\Delta t\mathbf{D}_{j}^{h}\in\mathbb{R}^{n_{x}\times n_{x}},$
$\displaystyle\mathbf{T}_{j}^{\text{imp}}$
$\displaystyle:=\left(\mathbf{I}-\Delta
t\mathbf{D}_{j}^{v}\right)^{-1}\in\mathbb{R}^{n_{x}\times n_{x}}$

for each time instant $t_{j}$, $j\in\\{0,\ldots,n_{t}-1\\}$, we obtain for a
time step of the marine ecosystem model using the TMM

$\displaystyle\mathbf{y}_{j+1}$
$\displaystyle=\mathbf{T}_{j}^{\text{imp}}\left(\mathbf{T}_{j}^{\text{exp}}\mathbf{y}_{j}+\Delta
t\mathbf{q}_{j}\left(\mathbf{y}_{j},\mathbf{u}\right)\right)=:\varphi_{j}\left(\mathbf{y}_{j},\mathbf{u}\right),$
$\displaystyle j$ $\displaystyle=0,\ldots,n_{t}-1.$ (8)

The transport matrices are sparse representations of the monthly averaged
tracer transport. The explicit and implicit transport matrices are sparse
because these matrices are generated with a grid-point based ocean circulation
model and the implicit ones (i.e., the inverse of the discretization matrices)
contain only the vertical part of the diffusion. In practical computations,
monthly averaged matrices are stored and interpolated linearly to compute an
approximation for any time instant $t_{j}$, $j=0,\ldots,n_{t}-1$, because
storing the matrices for all time-steps in a year is practically impossible.
Assuming annually periodicity of the ocean circulation, the TMM approximates
the ocean circulation using twelve pairs of pre-computed transport matrices.
In the present paper, these matrices are computed with the MIT ocean model
[31] using a global configuration with a latitudinal and longitudinal
resolution of $None$ and $15$ vertical layers.

###  2.4 Computation of steady annual cycles

A periodic solution (i.e., a steady annual cycle) in a fully discrete setting
fulfills

$\displaystyle\mathbf{y}_{n_{t}}$ $\displaystyle=\mathbf{y}_{0}$

applying the above iteration (8) over one year model time. The nonlinear
mapping

$\displaystyle\Phi$
$\displaystyle:=\varphi_{n_{t}-1}\circ\ldots\circ\varphi_{0}$

with $\varphi_{j}$ specified in (8) defines the time integration of (8) over
one model year. For a marine ecosystem model, the steady annual cycle, in
particular, is a fixed-point of this mapping. Starting from an arbitrary
vector $\mathbf{y}^{0}\in\mathbb{R}^{n_{y}n_{x}}$ and model parameters
$\mathbf{u}\in\mathbb{R}^{n_{u}}$, a classical fixed-point iteration takes the
form

$\displaystyle\mathbf{y}^{\ell+1}$
$\displaystyle=\Phi\left(\mathbf{y}^{\ell},\mathbf{u}\right),$
$\displaystyle\ell=0,1,\ldots.$ (9)

The vector $\mathbf{y}^{\ell}\in\mathbb{R}^{n_{y}n_{x}}$ contains the tracer
concentrations at the first time instant of the model year $\ell\in\mathbb{N}$
if we interpret the fixed-point iteration (9) as pseudo-time stepping or
_spin-up_ , a term which is widely used in ocean and climate research.

As high-fidelity model, we used $n_{f,\ell}:=3000$ model years to compute the
steady annual cycle with the spin-up. The accuracy of the spin-up is already
satisfactory after $3000$ model years [45, 23]. Nevertheless, the spin-up over
more than $3000$ model years can improve the accuracy. The vector

$\displaystyle\mathbf{y}_{f}\approx\left(\mathbf{y}^{\ell}_{k}\right)_{k=0,\ldots,n_{t}-1}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$

contains the trajectory of the steady annual cycle computed with the high-
fidelity model and includes the different tracer concentrations
$\mathbf{y}^{\ell}_{k}:=\varphi_{k-1}\circ\ldots\circ\varphi_{0}\left(\mathbf{y}^{\ell},\mathbf{u}\right)$,
for $k\in\\{1,\ldots,n_{t}-1\\}$, and
$\mathbf{y}^{\ell}_{0}:=\mathbf{y}^{\ell}$ for each time step in the model
year.

##  3 The artificial neural network

We applied the prediction of an ANN as an approximation of the steady annual
cycle of the N model [40]. For this purpose, we used a multilayer perceptron
with a feed-forward architecture consisting of an input layer, three hidden
layers and an output layer. More specifically, we employed a fully connected
network, i.e., there do exist all possible connections between the neurons in
two consecutive layers.

The ANN predicted the tracer concentration for the first time of the model
year using the $n_{u}=5$ model parameters as input. Therefore, the input layer
consisted of $5$ neurons, while the output layer contained $52\,749$ neurons.
A prediction of the steady annual cycle taking the whole trajectory into
account was not possible because the output layer would consist of almost
$2.4$ million neurons so that the training of the ANN would be infeasible. The
fully connected network consisted of three hidden layers with $10$, $25$ and
$211$ neurons, respectively.

We trained the ANN with supervised learning [17] using backpropagation and
used the stochastic gradient descent [48, 6] for the optimization of the
network weights. The loss function was the mean squared error. For the neurons
of the hidden layers, we applied the exponential linear unit [9] as activation
function and the scaled exponential linear unit [20] for the neurons of the
output layer.

The total mass of the predicted tracer concentration was essential for a
reasonable prediction because the marine ecosystem model preserves mass. We
adjusted the prediction $\mathbf{y}_{\text{ANN}}\in\mathbb{R}^{n_{x}}$ of the
ANN pointwise by

$\displaystyle\tilde{\mathbf{y}}_{\text{ANN}}$
$\displaystyle=\frac{\sum_{k=1}^{n_{x}}\left|V_{k}\right|m}{\sum_{k=1}^{n_{x}}\left|V_{k}\right|\mathbf{y}_{\text{ANN}}(x_{k},t_{0})}\,\mathbf{y}_{\text{ANN}}$

using the box volumes $\left|V_{k}\right|$, $k=1,\ldots,n_{x}$, of the
discretization and a global mean tracer concentration $m\in\mathbb{R}_{>0}$ in
order to obtain the required overall mass.

A given set of pairs of input and output data is divided into _training_ ,
_validation_ and _test data_ for the training procedure (i.e., the
optimization of the internal weights of the ANN) using supervised learning.
The optimization adjusts the internal weights by the use of the training data
and monitors the behavior of the already partially optimized network on the
validation data. The test data are used to assess the quality of the trained
network after finishing the training procedure.

##  4 Surrogate-based optimization

The SBO [3, 46, 15, 29] replaces the high- fidelity model by a surrogate
during the optimization. For many nonlinear optimization problems, the high
computational effort to evaluate the objective function is a major bottleneck.
Conversely, the SBO uses the surrogate to reduce the computational effort
because the evaluation of the surrogate is computationally much cheaper and
still represents a reasonable approximation of the high-fidelity model. In
order to obtain optimal model parameters close to the high-fidelity optimum,
the SBO keeps the surrogate close to the high-fidelity model by iterative
updates and re-optimization of the surrogate. As a result, the calculated
parameters of a well performing SBO are close to those of a conventional
optimization algorithm.

###  4.1 Cost function

For the optimization, we defined the discrete cost function

$\displaystyle J(\mathbf{z})$
$\displaystyle:=\frac{1}{2}\left\|\mathbf{z}-\mathbf{y}_{d}\right\|_{2}^{2}:=\frac{1}{2}\sum_{j=0}^{n_{t}-1}\sum_{i=1}^{n_{y}}\sum_{k=1}^{n_{x}}\left(z_{jik}-y_{d,jik}\right)^{2}$
(10)

by measuring the difference between an annual cycle
$\mathbf{z}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$, for example calculated during the
optimization process, and the target annual cycle
$\mathbf{y}_{d}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$ in an Euclidean norm. Both
vectors $\mathbf{z}$ and $\mathbf{y}_{d}$ are indexed as
$\mathbf{z}=\left(\left(\left(z_{jik}\right)_{k=1}^{n_{x}}\right)_{i=1}^{n_{y}}\right)_{j=0}^{n_{t}-1}$
and
$\mathbf{y}_{d}=\left(\left(\left(y_{jik}\right)_{k=1}^{n_{x}}\right)_{i=1}^{n_{y}}\right)_{j=0}^{n_{t}-1}$.
We did not use weights in the Euclidean norm because we performed the
following optimizations exclusively with synthetic data. When using real
measurement data, the cost function (10) can be restricted by weights, for
example, to take only discrete points in space and time for which measurements
exist or the error variances of the measurements into account. In addition,
the cost function can be restricted to individual tracers. [23] investigated
different cost functions for biogeochemical models.

###  4.2 Optimization

We used the SBO instead of a conventional optimization approach (for example
gradient-based or meta-heuristics) to identify the optimal parameters of a
biogeochemical model. These parameters are usually identified using a
nonlinear optimization problem

$\displaystyle\min_{\mathbf{u}\in
U_{ad}}J\left(\mathbf{y}_{f}(\mathbf{u})\right)$

with the set of admissible parameters
$U_{ad}:=\left\\{\mathbf{u}\in\mathbb{R}^{n_{u}}\leavevmode\nobreak\
:\leavevmode\nobreak\
\mathbf{b}_{\ell}\leq\mathbf{u}\leq\mathbf{b}_{u}\right\\}$ and parameter
boundaries $\mathbf{b}_{\ell},\mathbf{b}_{u}\in\mathbb{R}^{n_{u}}$ with
$\mathbf{b}_{\ell}\leq\mathbf{b}_{u}$. The inequalities with the parameter
boundaries are meant component-wise. We did not optimize the expensive high-
fidelity model ($\mathbf{y}_{f}$) using any conventional approach because the
computational effort is unaffordable.

The SBO replaces the high-fidelity model with its surrogate, whose evaluation
is computationally cheaper than the high-fidelity model but still provides a
reasonably accurate approximation of the high-fidelity model. There are
several ways of constructing the surrogate, such as _function-approximation_
surrogates [46, 51, 52] or _physically-based_ surrogates [53]. In the present
paper, we used two different methods to construct the low-fidelity model
($\mathbf{y}_{c}$): a physics-based low-fidelity model and a neural network.
These low-fidelity models are a less accurate but computationally more
efficient approximation of the high-fidelity model. However, the accuracy of
the low-fidelity models is not sufficient for a direct optimization. In order
to reduce the misalignment with respect to the high-fidelity model, we form
the surrogate from the low-fidelity model with a subsequent correction. Since
the surrogate approximates only the high-fidelity model adequately in a local
environment, we also used a trust-region safeguard [10, 21, 45] to limit the
step size in the optimization to a certain trusted region ensuring the
computation of a local minimum with the surrogate if both the low- and the
high-fidelity model are sufficiently smooth. This enhances any surrogate
optimization.

Data: $\mathbf{u}_{0}\in\mathbb{R}^{n_{u}}$, $\delta_{0}\in\mathbb{R}_{>0}$,
$\delta^{\text{min}}\in\mathbb{R}_{>0}$, $\gamma\in\mathbb{R}_{>0}$,
$r_{\text{decr}}\in\mathbb{R}_{>0}$, $r_{\text{incr}}\in\mathbb{R}_{>0}$,
$m_{\text{decr}}\in\mathbb{R}_{>0}$, $m_{\text{incr}}\in\mathbb{R}_{>0}$,
$a_{\ell}\in\mathbb{R}_{>0}$, $a_{u}\in\mathbb{R}_{>0}$,
$\delta\in\mathbb{R}_{>0}$,

Result: $\mathbf{u}_{k_{\text{min}}}\in\mathbb{R}^{n_{u}}$ with
$k_{\text{min}}=\min\left\\{k\in\mathbb{N}\leavevmode\nobreak\
:\leavevmode\nobreak\
\left\|\mathbf{u}_{k}-\mathbf{u}_{k-1}\right\|_{2}^{2}\leq\gamma\vee\delta_{k}\leq\delta_{k}^{\text{min}}\right\\}$

$k=0$

while
_$k=0\vee\left\|\mathbf{u}_{k}-\mathbf{u}_{k-1}\right\|_{2}^{2}\leq\gamma\vee\delta_{k}\leq\delta^{\text{min}}$_
do

Evaluate high-fidelity $\mathbf{y}_{f}(\mathbf{u}_{k})$ and low-fidelity model
$\mathbf{y}_{c}(\mathbf{u}_{k})$

// Construct surrogate

Compute correction vector $\mathbf{a}_{k}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$ with
$a_{k,jil}:=\begin{cases}1.0,&\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)_{jil}<\delta\vee\left(\mathbf{y}_{c}(\mathbf{u}_{k})\right)_{jil}<\delta\\\
a_{\ell},&\frac{\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)_{jil}}{\left(\mathbf{y}_{c}(\mathbf{u}_{k})\right)_{jil}}<a_{\ell}\\\
a_{u},&\frac{\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)_{jil}}{\left(\mathbf{y}_{c}(\mathbf{u}_{k})\right)_{jil}}>a_{u}\\\
\frac{\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)_{jil}}{\left(\mathbf{y}_{c}(\mathbf{u}_{k})\right)_{jil}},&\text{else}\end{cases}$,
$i\in\\{1,\ldots,n_{y}\\},j\in\\{0,\ldots,n_{t}-1\\},l\in\\{1,\ldots,n_{x}\\}$

$\mathbf{s}_{k}(\mathbf{u}):=\mathbf{a}_{k}\mathbf{y}_{c}(\mathbf{u})$

// Optimize surrogate

$\mathbf{u}_{k+1}=\operatorname*{argmin}_{\mathbf{u}\in
U_{ab},\left\|\mathbf{u}-\mathbf{u}_{k}\right\|_{2}\leq\delta_{k}}J(\mathbf{s}_{k}(\mathbf{u}))$

// Update trust-region radius

if  _$J\left(\mathbf{y}_{f}(\mathbf{u}_{k+1})\right)
<J\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)$ _ then

$\rho_{k}=\frac{J\left(\mathbf{y}_{f}(\mathbf{u}_{k+1})\right)-J\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)}{J\left(\mathbf{s}_{k}(\mathbf{u}_{k+1})\right)-J\left(\mathbf{s}_{k}(\mathbf{u}_{k})\right)}$

$\delta_{k+1}=\begin{cases}\frac{\delta_{k}}{m_{\text{decr}}},&\rho_{k}<r_{\text{decr}},\\\
\delta_{k},&r_{\text{decr}}\leq\rho_{k}\leq r_{\text{incr}}\\\ \delta_{k}\cdot
m_{\text{incr}},&\rho_{k}>r_{\text{incr}}\end{cases}$

$k=k+1$

else

$\delta_{k}=\frac{\delta_{k}}{m_{\text{decr}}}$

Algorithm 1 Surrogate-based optimization

The SBO consists of an iterative approach generating the surrogate again after
a successful surrogate optimization. The algorithm of the SBO shown in
Algorithm 1 can be divided into three parts. The first part contains the
evaluation of the low- and high-fidelity model including the construction of
the surrogate. The different low-fidelity models used in this paper are
explained in Section 4.2.1, 4.2.2 as well as 4.2.3, and the construction of
the surrogate as a correction of the low-fidelity model is presented in
Section 4.2.4. In the second part of the algorithm, the surrogate
($\mathbf{s}_{k}$) is optimized using a gradient-based algorithm, i.e.,

$\displaystyle\mathbf{u}_{k+1}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{u}\in
U_{ad},\left\|\mathbf{u}-\mathbf{u}_{k}\right\|\leq\delta_{k}}J\left(\mathbf{s}_{k}(\mathbf{u})\right)$

for the iteration $k\in\mathbb{N}_{0}$ and trust-region radius
$\delta_{k}\in\mathbb{R}_{>0}$. If the surrogate optimization reduces the cost
function value evaluated for the high-fidelity model, the parameter vector
$\mathbf{u}_{k+1}$ is accepted and the trust-region radius is updated
according to the rules below. Otherwise, the parameter vector
$\mathbf{u}_{k+1}$ is rejected, and the surrogate optimization is performed
again with the same surrogate but reduced trust-region radius. In analogy to
[45], we used the classical rules for updating the trust-region radius [10,
21] with slightly modified parameters, i.e.,

$\displaystyle\delta_{k+1}$
$\displaystyle=\begin{cases}\frac{\delta_{k}}{m_{\text{decr}}},&\rho_{k}<r_{\text{decr}}\\\
\delta_{k},&r_{\text{decr}}\leq\rho_{k}\leq r_{\text{incr}}\\\ \delta_{k}\cdot
m_{\text{incr}},&\rho_{k}>r_{\text{incr}}\end{cases}$

with the parameters $\delta_{0}=0.06$, $r_{\text{incr}}=0.75$,
$r_{\text{decr}}=0.01$ and $m_{\text{incr}},m_{\text{decr}}\in\mathbb{R}_{>0}$
as well as the gain ratio

$\displaystyle\rho_{k}$
$\displaystyle:=\frac{J\left(\mathbf{y}_{f}(\mathbf{u}_{k+1})\right)-J\left(\mathbf{y}_{f}(\mathbf{u}_{k})\right)}{J\left(\mathbf{s}_{k}(\mathbf{u}_{k+1})\right)-J\left(\mathbf{s}_{k}(\mathbf{u}_{k})\right)}.$

Consequently, the trust-region radius is decreased, on the one hand, if the
parameter vector $\mathbf{u}_{k+1}$ is rejected, or if the improvement of the
cost function value which was evaluated with the high-fidelity model is too
small compared to the prediction by the surrogate and increased, on the other
hand, if the prediction by the surrogate was adequate. The termination
condition of the iterative approach represents the third part of the
algorithm. We used a threshold $\gamma\in\mathbb{R}_{>0}$ for the absolute
step size in combination with a threshold
$\delta^{\text{min}}\in\mathbb{R}_{>0}$ for the trust-region radius.
Alternative termination conditions, for example, are certain convergence
criteria, an expected cost function value or a specific number of iterations.

####  4.2.1 Low-fidelity model using truncated spin-up

Our first low-fidelity model used the approach of a truncated spin-up (denoted
as $\mathbf{y}_{c}^{\text{N-DOP}}$ for the N-DOP model and
$\mathbf{y}_{c}^{\text{N}}$ for the N model) [45]. In order to compute an
approximation of the steady annual cycle using the low-fidelity model, we
reduced the number of model years to $n_{c}=25$ (cf. Section 2.4). The time
steps were, however, identical for the low- and high-fidelity model. The
number of model years $n_{f}=3000$ employed for the spin-up calculation of the
high-fidelity model is much larger than $n_{c}$ for the low-fidelity model.
Therefore, the computational effort for the low-fidelity model is
significantly lower than for the high-fidelity model.

####  4.2.2 Low-fidelity model using prediction of an ANN

Our second low-fidelity model approximated the steady annual cycle with the
prediction of an ANN [40]. For this purpose, the ANN predicted the tracer
concentration of a steady annual cycle approximation for the first time of the
model year. To receive the whole trajectory as low-fidelity model (denoted as
$\mathbf{y}_{c}^{\text{ANN}}$), we calculated the spin-up over one model year
using the prediction as initial concentration. We used the ANN described in
Section 3 including the adjustment of the prediction to obtain the required
overall mass. Particularly, the computational effort to evaluate this low-
fidelity model is very low because the computation of the trajectory for one
model year is already the main effort.

####  4.2.3 Low-fidelity model using ANN as initial value generator

For our third low-fidelity model, the prediction of the ANN served as initial
concentration for a truncated spin-up [40]. Instead of using the prediction
directly as low-fidelity model as in Section 4.2.2, we used the ANN described
in Section 3 to generate an initial concentration for the spin-up calculation
with a reduced number of model years $n_{c}=50$. The time step was the same as
for the high-fidelity model. In summary, the low-fidelity (denoted as
$\mathbf{y}_{c}^{\text{ANN-N}}$) model consisted of the ANN as initial value
generator in combination with a truncated spin-up. This low-fidelity model
reduced again the computationally effort with respect to the high-fidelity
model.

####  4.2.4 Multiplicative response correction

For the three low-fidelity models, the multiplicative response correction
reduces the misalignment of the low-fidelity model with respect to the high-
fidelity model. [43, 45] have already used this multiplicative correction
motivated by the physics-based approach of the low-fidelity model because the
steady annual cycle calculated with the low-fidelity model resembles that of
the high-fidelity model. For each iteration $k\in\mathbb{N}_{0}$, the
multiplicative correction vector
$\mathbf{a}_{k}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$ contains the pointwise ratio
of the high- and low-fidelity solution

$\displaystyle\mathbf{a}_{k}$
$\displaystyle:=\frac{\mathbf{y}_{f}(\mathbf{u}_{k})}{\mathbf{y}_{c}(\mathbf{u}_{k})}.$

In order to avoid larger entries in the correction vector, for example for
values of the low-fidelity model close to zero or for a value of the low-
fidelity model being several orders of magnitude smaller than the value of the
high-fidelity model, we restrict the entries of the correction vector by an
upper bound ($a_{u}\in\mathbb{R}_{>0}$) and a lower one
($a_{\ell}\in\mathbb{R}_{>0}$) and, furthermore, we set the value to $1.0$ if
the value of the low-fidelity model or high-fidelity model is close to zero,
i.e., falls below a threshold of $\delta=5\cdot 10^{-3}$. The pointwise
multiplication of the correction vector $\mathbf{a}_{k}$ and the low-fidelity
solution defines the surrogate, i.e.,

$\displaystyle\mathbf{s}_{k}\left(\mathbf{u}\right)$
$\displaystyle:=\mathbf{a}_{k}\mathbf{y}_{c}\left(\mathbf{u}\right)$

for $\mathbf{u}\in U_{ad}$. In particular, the construction of the surrogate
exclusively requires a single evaluation of the high-fidelity model. The
subsequent optimization needs only the evaluation of the surrogate and thus of
the low-fidelity model. Moreover, this construction of the surrogate improves
its accuracy and the performance of the entire SBO [45].

##  5 Results

We present in this section the results obtained with the three different low-
fidelity models. For the low-fidelity model described in Section 4.2.1, we
performed a parameter optimization for both biogeochemical models. In
contrast, we used the low-fidelity models described in Section 4.2.2 and 4.2.3
only for a parameter optimization using the N model because the neural network
in Section 3 was designed for the N model. In this paper, we have used
synthetic target data for the parameter optimization in order to compare the
different low-fidelity models. We assessed the optimization, on the one hand,
with the improvement of the cost function $J\left(\mathbf{y}_{f}\right)$ and,
on the other hand, using the accuracy of matching both the target phosphate
concentration and the optimal parameters.

###  5.1 Experimental setup

We configured the surrogate-based optimization uniformly for each low-fidelity
model. Although Table 2 lists the parameter vectors exclusively for the N-DOP
model, we used these parameter vectors also for the N model. For this purpose,
the parameter vector
$\left(k_{w},\mu_{P},K_{N},K_{I},\sigma_{\text{DOP}},\lambda^{\prime}_{\text{DOP}},b\right)$
of the N-DOP model had to be restricted to the parameter vector
$\left(k_{w},\mu_{P},K_{N},K_{I},b\right)$ of the N model, i.e., the
parameters $\sigma_{\text{DOP}}$ and $\lambda^{\prime}_{\text{DOP}}$ were
omitted.

Table 2: Parameter vectors for the biogeochemical models. Initial ($\mathbf{u}_{0}$) and optimal ($\mathbf{u}_{d}$) parameter vector used for the SBO runs, test parameter vector ($\bar{\mathbf{u}}$ and $\tilde{\mathbf{u}}$) to assess the quality of the surrogate model as well as lower $(\mathbf{b}_{\ell})$ and upper $(\mathbf{b}_{u})$ bounds for the parameter vectors. Parameter | $k_{w}$ | $\mu_{P}$ | $K_{N}$ | $K_{I}$ | $\sigma_{\text{DOP}}$ | $\lambda^{\prime}_{\text{DOP}}$ | $b$
---|---|---|---|---|---|---|---
$\mathbf{u}_{0}$ | 0.04 | 3.5 | 0.8 | 25.0 | 0.4 | 0.3 | 0.78
$\mathbf{u}_{d}$ | 0.02 | 2.0 | 0.5 | 30.0 | 0.67 | 0.5 | 0.858
$\bar{\mathbf{u}}$ | 0.016 | 1.6 | 0.4 | 24.0 | 0.536 | 0.4 | 0.686
$\tilde{\mathbf{u}}$ | 0.0193 | 1.93 | 0.483 | 29.0 | 0.648 | 0.483 | 0.829
$\mathbf{b}_{\ell}$ | 0.01 | 1.0 | 0.25 | 15.0 | 0.05 | 0.25 | 0.7
$\mathbf{b}_{u}$ | 0.05 | 4.0 | 1.0 | 60.0 | 0.95 | 1.0 | 1.5

####  5.1.1 Training of the ANN

We trained the fully connected network described in Section 3 using the
parameter vectors of a Latin hypercube sample [[, cf.]]McBeCo79. We created
the Latin hypercube sample with $1100$ parameter vectors within the bounds
($\mathbf{b}_{\ell}$ and $\mathbf{b}_{u}$) of Table 2 using the routine Lhs of
[28] and computed a steady annual cycle for each parameter vector. For the
computation of the steady annual cycle, we applied the marine ecosystem
toolkit for optimization and simulation in 3D (Metos3D) [42] and performed a
spin-up over $10\,000$ model years, starting from a global mean concentration
of $2.17\text{\,}\mathrm{mmol}\text{\,}\mathrm{P}\text{\,}{\mathrm{m}}^{-3}$
for $\textrm{PO}_{4}$. For the training of the fully connected network, we
used the parameter vectors of the Latin hypercube sample as input and the
tracer concentration of the first time instant in the last, the $10\,000$th,
model year as output. We trained the neural network over $1000$
iterations/epochs.

####  5.1.2 Setup of the surrogate-based optimization

For the calculation of a steady annual cycle, we applied Metos3D using
$n_{t}=45$ time steps. Except for the two low-fidelity models using the ANN
(see Section 4.2.2 and 4.2.3), the spin-up calculation started always with a
constant global mean concentration of
$2.17\text{\,}\mathrm{mmol}\text{\,}\mathrm{P}\text{\,}{\mathrm{m}}^{-3}$ for
$\textrm{PO}_{4}$ and, if present,
$0.0001\text{\,}\mathrm{mmol}\text{\,}\mathrm{P}\text{\,}{\mathrm{m}}^{-3}$
for DOP. For the synthetic data, we selected a random parameter vector
$\mathbf{u}_{d}\in\mathbb{R}^{n_{u}}$ (see Table 2) and computed the steady
annual cycle $\mathbf{y}_{d}\in\mathbb{R}^{n_{t}n_{y}n_{x}}$ using the spin-up
over $10\,000$ model years with $n_{t}=2880$ time steps and restricted the
steady annual cycle to the corresponding time steps for $n_{t}=45$, i.e.,

$\displaystyle\mathbf{y}_{d}:=\left(\mathbf{y}_{i}^{10000}\left(\mathbf{u}_{d}\right)\right)_{i\in\\{0,64,\ldots,2816\\}}.$

The SBO used the L-BFGS-B algorithm as optimization algorithm for the
surrogate [8, 58] implemented in SciPy111Scientific computing tools for
Python, https://www.scipy.org [55]. We started the SBO exemplary from a random
parameter vector $\mathbf{u}_{0}\in\mathbb{R}^{n_{u}}$ shown in Table 2.
Moreover, we limited the admissible parameter vectors using the lower
$\mathbf{b}_{\ell}$ and upper bound $\mathbf{b}_{u}$ (Table 2). The SBO
stopped if the absolute step size fell below a threshold of $\gamma=5\cdot
10^{-4}$, or the trust-region radius undercut a threshold of
$\delta^{\text{min}}=5\cdot 10^{-5}$. Finally, we restricted the entries of
the multiplicative correction vector for each surrogate using the lower and
upper bounds $a_{\ell}=0.1$ and $a_{u}=5$.

###  5.2 Suitability of the multiplicative correction technique

(a)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{y}_{c}^{\text{N}}(\bar{\mathbf{u}})$
(b)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{y}_{c}^{\text{ANN}}(\bar{\mathbf{u}})$
(c)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{y}_{c}^{\text{ANN-N}}(\bar{\mathbf{u}})$
(d)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{s}_{0}^{\text{N}}(\bar{\mathbf{u}})$
(e)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{s}_{0}^{\text{ANN}}(\bar{\mathbf{u}})$
(f)
$\mathbf{y}_{f}(\bar{\mathbf{u}})-\mathbf{s}_{0}^{\text{ANN-N}}(\bar{\mathbf{u}})$

(g)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{y}_{c}^{\text{N}}(\tilde{\mathbf{u}})$
(h)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{y}_{c}^{\text{ANN}}(\tilde{\mathbf{u}})$
(i)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{y}_{c}^{\text{ANN-N}}(\tilde{\mathbf{u}})$

(j)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{s}_{0}^{\text{N}}(\tilde{\mathbf{u}})$
(k)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{s}_{0}^{\text{ANN}}(\tilde{\mathbf{u}})$

(l)
$\mathbf{y}_{f}(\tilde{\mathbf{u}})-\mathbf{s}_{0}^{\text{ANN-N}}(\tilde{\mathbf{u}})$

Figure 1: Difference between the high-fidelity model $\mathbf{y}_{f}$ and the
low-fidelity model $\mathbf{y}_{c}$ as well as the surrogate $\mathbf{s}_{0}$
for the three different low-fidelity models (truncated-spin-up
$\mathbf{y}_{c}^{\text{N}}$ (left, Section 4.2.1), ANN prediction
$\mathbf{y}_{c}^{\text{ANN}}$ (middle, Section 4.2.2) and ANN prediction as
initial concentration for the truncated spin-up
$\mathbf{y}_{c}^{\text{ANN-N}}$ (right, Section 4.2.3) as low-fidelity model).
For the two parameter vectors $\bar{\mathbf{u}}$ and $\tilde{\mathbf{u}}$, the
phosphate concentration difference are shown on the surface
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) at the first time step of the
model year (in January). The surrogate was built using the parameter vector
$\mathbf{u}_{d}$, respectively.

The multiplicative response correction was suitable for the three low-fidelity
models. Figure 1 indicates the improvement of the low-fidelity model accuracy
using the multiplicative response correction. Although the concentrations of
the low-fidelity models differed significantly from each other in some areas,
the approximations of the high-fidelity model were very similar after the
multiplicative correction. Particularly, the surrogate approximated the high-
fidelity model better for each low-fidelity model if the parameter vector was
closer to the parameter vector used to construct the surrogate. The
approximation was, nevertheless, slightly worse using the prediction of the
ANN as low-fidelity model while the other two low-fidelity models showed only
marginal differences in the approximations after the correction. Except for
the low-fidelity model using the ANN prediction as detailed in Section 5.4,
the surrogates provided an appropriate approximation of the high-fidelity
model in the neighborhood of the parameter vector used to construct the
surrogate.

###  5.3 Low-fidelity model as truncated spin-up

In this section, we present the parameter identification with the SBO using
the low-fidelity model based on the truncated spin-up (see Section 4.2.1) for
both biogeochemical models. For the update of the trust-region radius, we
applied the parameter $m_{\text{incr}}=3$ for both biogeochemical models, as
well as the parameter $m_{\text{decr}}=20$ for the N-DOP model and
$m_{\text{decr}}=3$ for the N model.

####  5.3.1 N-DOP model

(a) Cost function value $J\left(\mathbf{y}_{f}\right)$. (b) Step-size norm of
the trust-region radius $\delta_{k}$.

Figure 2: Convergence of the cost function value and the step-size norm of the
trust-region radius for the exemplary SBO run with the truncated spin-up as
low-fidelity model ($\mathbf{y}_{c}^{\text{N-DOP}}$) and the N-DOP model.

(a)  (b)  (c)  (d)  (e)  (f)  (g)

Figure 3: Convergence of the single parameter values $\mathbf{u}_{k,i}$ for
each iteration of the exemplary SBO run with the truncated spin-up as low-
fidelity model ($\mathbf{y}_{c}^{\text{N-DOP}}$) and the N-DOP model.

The parameter identification using the SBO identified six out of seven model
parameters of the N-DOP model. The steady annual cycles calculated for each
iteration of the SBO converged to the target concentration $\mathbf{y}_{d}$
(Figure 2). Furthermore, the single parameter values converged against the
values of the optimal parameter vector $\mathbf{u}_{d}$ (see Figure 3), except
for parameter $K_{I}$. Especially in the iterations four to eight, the SBO
reduced the difference to the optimal parameter vector whereas even the step
size increased significantly. The SBO terminated after ten iterations because
the step size fell below the threshold of $\gamma=5\cdot 10^{-5}$ (i.e.,
$\left\|\mathbf{u}_{10}-\mathbf{u}_{9}\right\|_{2}^{2}<\gamma$). Although the
main patterns of the target phosphate concentration are already visible for
the solution $\mathbf{y}_{f}(\mathbf{u}_{0})$ of the initial parameter vector
$\mathbf{u}_{0}$, Figures 4 and 5 show the further improvement of the match
(for instance in the North Atlantic on the surface layer) during the
optimization.

(a) $\mathbf{y}_{d}$ (b) $\mathbf{y}_{f}(\mathbf{u}_{0})$ (c)
$\mathbf{y}_{f}(\mathbf{u}_{2})$ (d) $\mathbf{y}_{f}(\mathbf{u}_{5})$ (e)
$\mathbf{y}_{f}(\mathbf{u}_{8})$ (f) $\mathbf{y}_{f}(\mathbf{u}_{10})$

Figure 4: High-fidelity model output $\mathbf{y}_{f}$ obtained by the
exemplary SBO run with the truncated spin-up as low-fidelity model
($\mathbf{y}_{c}^{\text{N-DOP}}$) and the N-DOP model at the beginning and
after two, five, eight and ten iterations (i.e., evaluated with the parameter
vectors $\mathbf{u}_{0}$, $\mathbf{u}_{2}$, $\mathbf{u}_{5}$, $\mathbf{u}_{8}$
and $\mathbf{u}_{10}$) as well as the target data $\mathbf{y}_{d}$. Shown are
the phosphate concentrations on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) at the first time step of the
model year (in January).

(a) Location $None$ W, $None$ N (b) Location $None$ E, $None$ N

Figure 5: Annual cycles of phosphate of the high-fidelity model output
$\mathbf{y}_{f}$ obtained by the exemplary SBO run with the truncated spin-up
as low-fidelity model ($\mathbf{y}_{c}^{\text{N-DOP}}$) and the N-DOP model at
the beginning and after two, five, eight and ten iterations (i.e., evaluated
with the parameter vectors $\mathbf{u}_{0}$, $\mathbf{u}_{2}$,
$\mathbf{u}_{5}$, $\mathbf{u}_{8}$ and $\mathbf{u}_{10}$) as well as the
target data $\mathbf{y}_{d}$ on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) for two distinct locations.

The SBO was computationally very efficient. The exemplary SBO run required
only $15$ high-fidelity model evaluations, whereby the surrogate optimization
of four iterations were discarded due to a larger cost function value after
the optimization. In contrast, the surrogate optimization needed altogether
more than $6000$ evaluations of the surrogate. This large number of surrogate
evaluations resulted from the limitation of the iterations of the surrogate
optimization to $100$. Limiting the number of iterations to $10$, the SBO run
converged with a cost function value of about $650$ (instead of about $100$
using a limit of $100$ iterations as illustrated by Figure 2a) but with a
significantly lower number of about $1500$ surrogate evaluations.

####  5.3.2 N model

(a) Cost function value $J\left(\mathbf{y}_{f}\right)$. (b) Step-size norm of
the trust-region radius $\delta_{k}$.

Figure 6: Convergence of the cost function value and the step-size norm of the
trust-region radius for the exemplary SBO run with the truncated spin-up as
low-fidelity model ($\mathbf{y}_{c}^{\text{N}}$) and the N model.

(a)  (b)  (c)  (d)  (e)

Figure 7: Convergence of the single parameter values $\mathbf{u}_{k,i}$ for
each iteration of the exemplary SBO run with the truncated spin-up as low-
fidelity model ($\mathbf{y}_{c}^{\text{N}}$) and the N model.

(a) $\mathbf{y}_{d}$ (b) $\mathbf{y}_{f}(\mathbf{u}_{0})$ (c)
$\mathbf{y}_{f}(\mathbf{u}_{4})$ (d) $\mathbf{y}_{f}(\mathbf{u}_{7})$ (e)
$\mathbf{y}_{f}(\mathbf{u}_{10})$ (f) $\mathbf{y}_{f}(\mathbf{u}_{13})$

Figure 8: High-fidelity model output $\mathbf{y}_{f}$ obtained by the
exemplary SBO run with the truncated spin-up as low-fidelity model
($\mathbf{y}_{c}^{\text{N}}$) and the N model at the beginning and after $4$,
$7$, $10$ and $13$ iterations (i.e., evaluated with the parameter vectors
$\mathbf{u}_{0}$, $\mathbf{u}_{4}$, $\mathbf{u}_{7}$, $\mathbf{u}_{10}$ and
$\mathbf{u}_{13}$) as well as the target data $\mathbf{y}_{d}$. Shown are the
phosphate concentrations on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) at the first time step of the
model year (in January).

(a) Location $None$ W, $None$ N (b) Location $None$ E, $None$ N

Figure 9: Annual cycles of phosphate of the high-fidelity model output
$\mathbf{y}_{f}$ obtained by the exemplary SBO run with the truncated spin-up
as low-fidelity model ($\mathbf{y}_{c}^{\text{N}}$) and the N model at the
beginning and after $4$, $7$, $10$ and $13$ iterations (i.e., evaluated with
the parameter vectors $\mathbf{u}_{0}$, $\mathbf{u}_{4}$, $\mathbf{u}_{7}$,
$\mathbf{u}_{10}$ and $\mathbf{u}_{13}$) as well as the target data
$\mathbf{y}_{d}$ on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) for two distinct locations.

The parameters converged to the optimal parameter $\mathbf{u}_{d}$ using the
SBO for the N model with the exception of parameter $K_{I}$. The cost function
values $J(\mathbf{y}_{f})$ in Figure 6a indicate a solution close to the
target data $\mathbf{y}_{d}$. Particularly, the SBO identified four out of
five model parameters as detailed in Figure 7. Although the first iteration
already determined the parameter $b$, the SBO required twelve iterations to
identify the other three parameters. Thereby, the algorithm repeatedly
increased the step size (Figure 6b). As a result of falling below the
threshold $\delta^{\text{min}}$ for the trust-region radius, the SBO
terminated after $13$ iterations. According to Figures 8 and 9, the main
patterns of the target phosphate concentration are apparent for the tracer
concentrations of the selected parameter vectors. However, the SBO improved
the match sequentially (for example in the Indian Ocean, see also Figure 6a).
Nevertheless, local worsening occurred occasionally compared to the previously
calculated concentration (cf. Figure 9).

The computational costs of the SBO were low. The exemplary SBO run required
only $34$ high-fidelity model evaluations and about $6000$ surrogate
evaluations. Compared to the small number of $13$ iterations, more than half
of the high-fidelity model evaluations resulted from the calculation of the
reference cost function value for parameter vectors that had been calculated
by the surrogate optimization but did not cause a reduction of the cost
function. Analogous to the SBO using the N-DOP model in the previous section
5.3.1, the number of surrogate evaluations could be reduced by limiting the
optimization iterations of the surrogate optimization to $10$. As a
consequence, the final cost function value increased significantly to over
$1000$ for such an optimization run compared to a value of about $100$ when a
limit of $100$ iterations of the surrogate optimization was used as shown in
Figure 6a.

###  5.4 Prediction of an ANN as low-fidelity model

(a) Location $None$ W, $None$ N (b) Location $None$ E, $None$ N

Figure 10: Annual cycle of phosphate for the high- ($\mathbf{y}_{f}$), low-
fidelity model ($\mathbf{y}_{c}^{\text{ANN}}$) and surrogate
($\mathbf{s}_{0}$) on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) for two distinct locations.
Shown are the annual cycles for three parameter values $\mathbf{u}_{0}$ (as
“reference point”), $\bar{\mathbf{u}}$ (a neighbouring point) and
$\tilde{\mathbf{u}}$ (a point in a closer vicinity). The surrogate was built
using parameter vector $\mathbf{u}_{0}$, the reason why the solution of the
surrogate is omitted at the reference point.

Using the prediction of the ANN defined in Section 3 as low-fidelity model
(see Section 4.2.2) was unsuitable to construct a reliable surrogate in
conjunction with the multiplicative response correction for the SBO with the N
model. Figure 10 demonstrates the nearly identical phosphate concentrations
predicted by the ANN for the three parameter vectors
$\mathbf{u}_{0},\bar{\mathbf{u}},\tilde{\mathbf{u}}\in\mathbb{R}^{n_{u}}$
defined in Table 2. Consequently, the steady annual cycles calculated as spin-
up over one model year using the prediction as an initial concentration
differed only marginally from each other due to the different model
parameters. Instead of the high-fidelity solution of the corresponding
parameter vector, the surrogate, therefore, in each case, approximated the
high-fidelity solution for parameter vector $\mathbf{u}_{0}$ with which
surrogate $\mathbf{s}_{0}$ was constructed. We observed a similar behavior
using random parameter vectors in a small area around the reference parameters
used to construct the surrogate.

Figure 11: Convergence of the cost function value
$J\left(\mathbf{y}_{f}\right)$ for the exemplary SBO run with the N model
using the prediction of the ANN as low-fidelity model
($\mathbf{y}_{c}^{\text{ANN}}$). The iterations $0$ and $1$ (with parameter
vectors $\mathbf{u}_{0}$ and $\mathbf{u}_{1}$) of the SBO run correspond to
$1$ and $4$ evaluations of the high-fidelity model, respectively.

(a) $\mathbf{y}_{d}$ (b) $\mathbf{y}_{f}(\mathbf{u}_{0})$ (c)
$\mathbf{y}_{f}(\mathbf{u}_{1})$

Figure 12: High-fidelity model output $\mathbf{y}_{f}$ obtained by the
exemplary SBO run with the N model using the prediction of the ANN as low-
fidelity model ($\mathbf{y}_{c}^{\text{ANN}}$) at the beginning and after one
iteration (i.e., evaluated with the parameter vectors $\mathbf{u}_{0}$ and
$\mathbf{u}_{1}$) as well as the target data $\mathbf{y}_{d}$. Shown are the
phosphate concentrations on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) at the first time step of the
model year (in January).

The SBO did not identify any model parameters. After three surrogate
optimizations, the SBO calculated with surrogate $\mathbf{s}_{0}$ just one
parameter vector $\mathbf{u}_{1}\in\mathbb{R}^{n_{u}}$ for which the
corresponding reference cost function value was smaller than for the initial
parameter vector $\mathbf{u}_{0}$ (Figure 11). The subsequent surrogate
optimizations with surrogate $\mathbf{s}_{1}$ did not yield a parameter vector
whose high-fidelity model solution approximated the target data better.
Although the main pattern of the target phosphate concentration for the high-
fidelity solution corresponding to the parameter vector $\mathbf{u}_{0}$ as
well as $\mathbf{u}_{1}$ were evident as seen in Figure 12, there were many
differences, such as in the Pacific, and, consequently, parameter vector
$\mathbf{u}_{1}$ did not provide a satisfactory identification of the optimal
parameter vector $\mathbf{u}_{d}$. After only one iteration and ten high-
fidelity model evaluations, the SBO terminated because the trust-region radius
was decreased in every step using $m_{\text{decr}}=4$ and, hence, fell below
the threshold $\delta^{\text{min}}$ after ten steps.

###  5.5 ANN as initial value generator for the low-fidelity model

(a) Generation of initial values with an ANN for a spin-up over $50$ model
years as low-fidelity model ($\mathbf{y}_{c}^{\text{ANN-N}}$). (b) Spin-up
over $25$ model years using constant initial values as low-fidelity model
($\mathbf{y}_{c}^{\text{N}}$).

Figure 13: Annual cycle of phosphate for the high- ($\mathbf{y}_{f}$), low-
fidelity model ($\mathbf{y}_{c}$) and surrogate ($\mathbf{s}_{0}$) on the
surface layer ($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) for the location
$None$ W, $None$ N. Shown are the annual cycles for three parameter values
$\mathbf{u}_{0}$ (as “reference point”), $\bar{\mathbf{u}}$ (a neighbouring
point) and $\tilde{\mathbf{u}}$ (a point in a closer vicinity) defined in
Table 2. The surrogate is built using parameter vector $\mathbf{u}_{0}$, the
reason why the solution of the surrogate is omitted at the reference point.

Using the prediction of the ANN as initial concentration for a spin-up over
$50$ model years as low-fidelity model (see Section 4.2.3) was reasonable to
construct a reliable surrogate in conjunction with the multiplicative response
correction for the SBO using the N model. The annual cycles of phosphate in
Figure 13 indicate an almost equivalent approximation of the high-fidelity
model concentration by the surrogates using, on the one hand, the low-fidelity
model with initial concentration predicted by the ANN and, on the other hand,
the low-fidelity model with constant initial concentrations. Although the
concentrations differed for both low-fidelity models, the multiplicative
correction technique improved the accuracy of the low-fidelity models and
resulted in nearly the same approximation.

The parameter identification using the SBO identified four out of five model
parameters of the N model. Figure 14 shows the convergence of the tracer
concentrations to the target data $\mathbf{y}_{d}$ using different step sizes
during the optimization. Except for parameter $K_{I}$, the single parameters
converged to the parameter values of the optimal parameter vector
$\mathbf{u}_{d}$ (Figure 15). As a result of a step size less than the
threshold $\gamma$, the SBO terminated after $20$ iterations. For this
optimization run, we applied parameters $m_{\text{incr}}=2$ and
$m_{\text{decr}}=2$ for the update of the trust-region radius. Figures 16 and
17 illustrate the continuous improvement of the agreement of the high-fidelity
solution and the target data during the optimization.

(a) Cost function value $J\left(\mathbf{y}_{f}\right)$. (b) Step-size norm of
the trust-region radius $\delta_{k}$.

Figure 14: Convergence of the cost function value and the step-size norm of
the trust-region radius for the exemplary SBO run with the prediction of the
ANN as initial concentration for the spin-up as low-fidelity model
($\mathbf{y}_{c}^{\text{ANN-N}}$).

(a)  (b)  (c)  (d)  (e)

Figure 15: Convergence of the single parameter values $\mathbf{u}_{k,i}$ for
each iteration of the exemplary SBO run with the prediction of the ANN as
initial concentration for the spin-up as low-fidelity model
($\mathbf{y}_{c}^{\text{ANN-N}}$).

(a) $\mathbf{y}_{d}$ (b) $\mathbf{y}_{f}(\mathbf{u}_{0})$ (c)
$\mathbf{y}_{f}(\mathbf{u}_{5})$ (d) $\mathbf{y}_{f}(\mathbf{u}_{10})$ (e)
$\mathbf{y}_{f}(\mathbf{u}_{15})$ (f) $\mathbf{y}_{f}(\mathbf{u}_{20})$

Figure 16: High-fidelity model output $\mathbf{y}_{f}$ obtained by the
exemplary SBO run with the prediction of the ANN as initial concentration for
the spin-up as low-fidelity model ($\mathbf{y}_{c}^{\text{ANN-N}}$) at the
beginning and after $5$, $10$, $15$ and $20$ iterations (i.e., evaluated with
the parameter vectors $\mathbf{u}_{0}$, $\mathbf{u}_{5}$, $\mathbf{u}_{10}$,
$\mathbf{u}_{15}$ and $\mathbf{u}_{20}$) as well as the target data
$\mathbf{y}_{d}$. Shown are the phosphate concentrations on the surface layer
($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) at the first time step of the
model year (in January).

(a) Location $None$ W, $None$ N (b) Location $None$ E, $None$ N

Figure 17: Annual cycles of phosphate of the high-fidelity model output
$\mathbf{y}_{f}$ obtained by the exemplary SBO run with the prediction of the
ANN as initial concentration for the spin-up as low-fidelity model
($\mathbf{y}_{c}^{\text{ANN-N}}$) at the beginning and after $5$, $10$, $15$
and $20$ iterations (i.e., evaluated with the parameter vectors
$\mathbf{u}_{0}$, $\mathbf{u}_{5}$, $\mathbf{u}_{10}$, $\mathbf{u}_{15}$ and
$\mathbf{u}_{20}$) as well as the target data $\mathbf{y}_{d}$ on the surface
layer ($0\text{\,}\mathrm{m}50\text{\,}\mathrm{m}$) for two distinct
locations.

The computational effort of the parameter identification using the SBO was
low. To obtain parameter vector $\mathbf{u}_{20}$, the SBO required $35$ high-
fidelity model evaluations and about $8000$ evaluations of the surrogate. In
the exemplary SBO run, the surrogate optimization needed up to $50$
iterations. Indeed, the number of surrogate evaluations could be reduced if we
restricted the iterations for each surrogate optimization to ten. Apart from
reducing computational costs, the limitation of the iterations did not lead to
an adequate identification of the model parameters because the cost function
reached a value of $1300$ only.

##  6 Conclusions

The SBO is based on a low-fidelity model that, on the one hand, can be
evaluated computationally inexpensive and, on the other hand, is an
appropriate approximation of the high-fidelity model. A suitable correction
technique improves additionally the accuracy of the low-fidelity model
building the surrogate. We have compared three different low-fidelity models
in conjunction with the multiplicative response correction already employed by
[45]. In addition to the truncated spin-up [45], we applied a neural network
only and in combination with the truncated spin-up as low-fidelity model.

Using the SBO, the parameter optimization yielded similar results for the low-
fidelity models including and excluding the application of the neural network
except for the direct use of the prediction. With the exception of model
parameter $K_{I}$ modeling the light intensity compensation, the SBO
identified all model parameters for both biogeochemical models and confirmed
the results of [45] for the N-DOP model. However, the parameter
identifications were not optimal. A reason for the mismatch, on the one hand,
was possibly the low accuracy of the underlying low-fidelity model as well as
a weak sensitivity with regard to the model parameters as already observed by
[47] for a different biogeochemical model. For instance, the construction of
the surrogate considering information about the high-fidelity model could
increase the accuracy. On the other hand, the trust region radius limited the
convergence of parameter $K_{I}$ because the admissible parameter range of
this parameter was much larger than for the other parameters and, thus, larger
adjustments were necessary.

The use of neural networks to construct the surrogate reduced significantly
the computational costs of a parameter optimization. A reduction of these
costs is desirable because a parameter identification in climate models is
very expensive [25, 26, 33]. Indeed, the SBO, however, reduces the
computational costs compared to standard optimization routines but the
reduction of costs for the low-fidelity model decreases these further. If a
neural network is available, the computational costs of the prediction are
negligible. Therefore, the costs only originated from the subsequent truncated
spin-up, if this was used. However, the prediction of the ANN alone was not
suitable as low-fidelity model because only minimal concentration changes
arose from the prediction using various parameter vectors in a small
environment. A neural network trained with parameter vectors in such a small
environment, nevertheless, had the ability to predict the concentrations
reasonably for a random parameter vector of this environment whereas the
prediction was often insufficient for parameter vectors outside this
environment. The training of a neural network, furthermore, is computationally
costly and depends on the availability of a sufficient amount of training
data. Future work should, therefore, include the training of an ANN that
better predicts the concentrations locally and, thus, serves as low-fidelity
model for the construction of a reliable surrogate with the multiplicative
correction. For this purpose, reinforcement learning [54, 36] could help to
continuously improve the neural network during an optimization run. Moreover,
further future work should include the application of the SBO for parameter
identification using a normalization of the admissible parameter range for
each component of the parameter vector as well as using real measurement data
that, for example, is necessary for model assessment [[, see e.g.,]]KrKhOs10
and model calibration [[, see e.g.,]]KSKSO17, Kri17.

In summary, the main points of this paper are the following:

  * •

Using the truncated spin-up as a low-fidelity model, the SBO was a
computationally efficient method for the parameter identification of marine
ecosystem models.

  * •

The prediction of the presented ANN as low-fidelity model was unsuitable for
an SBO.

  * •

Using the prediction of the ANN as an initial concentration for a truncated
spin-up, the low-fidelity model represented an acceptable alternative to the
truncated spin-up with slightly higher computational costs.

## Code and data availability

The code used to generate the data in this publication is available at
https://github.com/slawig/bgc-ann and https://metos3d.github.io/. All used and
generated data are available at https://doi.org/10.5281/zenodo.5643667 [41].

## References

  * [1] Robert B. Bacastow and Ernst Maier-Reimer  “Ocean-circulation model of the carbon cycle”  In _Climate Dynamics_ 4.2, 1990, pp. 95–125  DOI: 10.1007/BF00208905
  * [2] Robert B. Bacastow and Ernst Maier-Reimer  “Dissolved organic carbon in modeling oceanic new production”  In _Global Biogeochemical Cycles_ 5.1, 1991, pp. 71–85  DOI: 10.1029/91GB00015
  * [3] John W. Bandler et al.  “Space mapping: the state of the art”  In _IEEE Transactions on Microwave Theory and Techniques_ 52.1, 2004, pp. 337–361  DOI: 10.1109/TMTT.2003.820904
  * [4] Harvey Thomas Banks and Karl Kunisch  “Estimation Techniques for Distributed Parameter Systems”, Systems and Control: Foundations and Applications  Basel: Birkhäuser, 1989  DOI: 10.1007/978-1-4612-3700-6
  * [5] Erik Bernsen, Henk A. Dijkstra and Fred W. Wubs  “A method to reduce the spin-up time of ocean models”  In _Ocean Modelling_ 20.4, 2008, pp. 380–392  DOI: 10.1016/j.ocemod.2007.10.008
  * [6] Léon Bottou and Olivier Bousquet  “The Tradeoffs of Large Scale Learning”  In _Advances in Neural Information Processing Systems_ 20, NIPS’07  Curran Associates Inc., 2008, pp. 161–168
  * [7] Kirk Bryan  “Accelerating the Convergence to Equilibrium of Ocean-Climate Models”  In _Journal of Physical Oceanography_ 14.4, 1984, pp. 666–673  DOI: 10.1175/1520-0485(1984)014<0666:ATCTEO>2.0.CO;2
  * [8] Richard H. Byrd, Peihuang Lu, Jorge Nocedal and Ciyou Zhu  “A Limited Memory Algorithm for Bound Constrained Optimization”  In _SIAM Journal on Scientific Computing_ 16.5, 1995, pp. 1190–1208  DOI: 10.1137/0916069
  * [9] Djork-Arné Clevert, Thomas Unterthiner and Sepp Hochreiter  “Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs)”, Preprint on arXiv, 2015  arXiv:1511.07289 [cs.LG]
  * [10] Andrew R. Conn, Nicholas I.. Gould and Philippe L. Toint  “Trust Region Methods”, MOS-SIAM Series on Optimization  Society for IndustrialApplied Mathematics, 2000  DOI: 10.1137/1.9780898719857
  * [11] Gokhan Danabasoglu, James C. McWilliams and William G. Large  “Approach to Equilibrium in Accelerated Global Oceanic Models”  In _Journal of Climate_ 9.5, 1996, pp. 1092–1110  DOI: 10.1175/1520-0442(1996)009<1092:ATEIAG>2.0.CO;2
  * [12] “Ocean Biogeochemistry”, Global Change – The IGBP Series  Berlin et al.: Springer, 2003  DOI: 10.1007/978-3-642-55844-3
  * [13] Katja Fennel, Martin Losch, Jens Schröter and Manfred Wenzel  “Testing a marine ecosystem model: sensitivity analysis and parameter optimization”  In _Journal of Marine Systems_ 28.1, 2001, pp. 45–63  DOI: 10.1016/S0924-7963(00)00083-X
  * [14] Wolfgang Fennel and Thomas Neumann  “Introduction to the Modelling of Marine Ecosystems” 72, Elsevier Oceanography Series  Elsevier, 2004
  * [15] Alexander I.J. Forrester and Andy J. Keane  “Recent advances in surrogate-based optimization”  In _Progress in Aerospace Sciences_ 45.1, 2009, pp. 50–79  DOI: 10.1016/j.paerosci.2008.11.001
  * [16] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning”, Adaptive computation and machine learning  MIT Press, 2016
  * [17] Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The elements of statistical learning”, Springer Series in Statistics  New York, NY: Springer, 2009, pp. 745  DOI: 10.1007/978-0-387-84858-7
  * [18] Samar Khatiwala  “A computational framework for simulation of biogeochemical tracers in the ocean”  In _Global Biogeochemical Cycles_ 21.3, 2007  DOI: 10.1029/2007GB002923
  * [19] Samar Khatiwala, Martin Visbeck and Mark A. Cane  “Accelerated simulation of passive tracers in ocean circulation models”  In _Ocean Modelling_ 9.1 Elsevier, 2005, pp. 51–69  DOI: 10.1016/j.ocemod.2004.04.002
  * [20] Günter Klambauer, Thomas Unterthiner, Andreas Mayr and Sepp Hochreiter  “Self-Normalizing Neural Networks”  In _Advances in Neural Information Processing Systems_ 30 Curran Associates, Inc., 2017, pp. 971–980
  * [21] Slawomir Koziel, John W. Bandler and Qingsha S. Cheng  “Robust Trust-Region Space-Mapping Algorithms for Microwave Design Optimization”  In _IEEE Transactions on Microwave Theory and Techniques_ 58.8, 2010, pp. 2166–2174  DOI: 10.1109/TMTT.2010.2052666
  * [22] Iris Kriest  “Calibration of a simple and a complex model of global marine biogeochemistry”  In _Biogeosciences_ 14.21, 2017, pp. 4965–4984  DOI: 10.5194/bg-14-4965-2017
  * [23] Iris Kriest, Samar Khatiwala and Andreas Oschlies  “Towards an assessment of simple global marine biogeochemical models of different complexity”  In _Progress In Oceanography_ 86.3-4, 2010, pp. 337–360  DOI: 10.1016/j.pocean.2010.05.002
  * [24] Iris Kriest et al.  “Calibrating a global three-dimensional biogeochemical ocean model (MOPS-1.0)”  In _Geoscientific Model Development_ 10.1, 2017, pp. 127–154  DOI: 10.5194/gmd-10-127-2017
  * [25] Eun Young Kwon and François Primeau  “Optimization and sensitivity study of a biogeochemistry ocean model using an implicit solver and in situ phosphate data”  In _Global Biogeochemical Cycles_ 20.4, 2006  DOI: 10.1029/2005GB002631
  * [26] Eun Young Kwon and François Primeau  “Optimization and sensitivity of a global biogeochemistry ocean model using combined in situ DIC, alkalinity, and phosphate data”  In _Journal of Geophysical Research: Oceans_ 113.C8, 2008  DOI: 10.1029/2007JC004520
  * [27] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _Nature_ 521, 2015, pp. 436–444  DOI: 10.1038/nature14539
  * [28] Abraham Lee  “pyDOE: Design of Experiments for Python” available at: https://pythonhosted.org/pyDOE/index.html (last access: November 30, 2021), 2014
  * [29] Leifur Leifsson and Slawomir Koziel  “Multi-fidelity design optimization of transonic airfoils using shape-preserving response prediction” ICCS 2010  In _Procedia Computer Science_ 1.1, 2010, pp. 1311–1320  DOI: 10.1016/j.procs.2010.04.146
  * [30] Alfred James Lotka  “Contribution to the Theory of Periodic Reactions”  In _The Journal of Physical Chemistry_ 14.3, 1910, pp. 271–274  DOI: 10.1021/j150111a004
  * [31] John Marshall et al.  “A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers”  In _Journal of Geophysical Research: Oceans_ 102.C3, 1997, pp. 5753–5766  DOI: 10.1029/96JC02775
  * [32] John H. Martin, George A. Knauer, David M. Karl and William W. Broenkow  “VERTEX: carbon cycling in the northeast Pacific”  In _Deep Sea Research Part A. Oceanographic Research Papers_ 34.2, 1987, pp. 267–285  DOI: 10.1016/0198-0149(87)90086-0
  * [33] Jann Paul Mattern, Katja Fennel and Michael Dowd  “Estimating time-dependent parameters for a biological ocean model using an emulator approach”  In _Journal of Marine Systems_ 96-97, 2012, pp. 32–47  DOI: 10.1016/j.jmarsys.2012.01.015
  * [34] Kendal McGuffie and Ann Henderson-Sellers  “The Climate Modelling Primer”  Chichester, West Sussex, United Kingdom: John Wiley & Sons, 2014
  * [35] Michael D. McKay, Richard J. Beckman and William Jay Conover  “A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code”  In _Technometrics_ 21.2 JSTOR, 1979, pp. 239–245  DOI: 10.2307/1268522
  * [36] Volodymyr Mnih et al.  “Human-level control through deep reinforcement learning”  In _Nature_ 518.7540, 2015, pp. 529–533  DOI: 10.1038/nature14236
  * [37] Andreas Oschlies  “On the use of data assimilation in biogeochemical modelling”  In _Ocean Weather Forecasting_ Dordrecht: Springer, 2006, pp. 525–547  DOI: 10.1007/1-4020-4028-8_24
  * [38] Garth William Paltridge and C… Platt  “Radiative Processes in Meteorology and Climatology”  New York: Elsevier, 1976  DOI: 10.1002/qj.49710343713
  * [39] P. Parekh, Michael J. Follows and Edward A. Boyle  “Decoupling of iron and phosphate in the global ocean”  In _Global Biogeochemical Cycles_ 19.2, 2005  DOI: 10.1029/2004GB002280
  * [40] Markus Pfeil and Thomas Slawig  “Approximation of a marine ecosystem model by artificial neural networks”  In _Electronic Transactions on Numerical Analysis_ , 2021
  * [41] Markus Pfeil and Thomas Slawig  “Surrogate-based optimization using an artificial neural network for a parameter identification in a 3D marine ecosystem model”  Zenodo, 2021  DOI: 10.5281/zenodo.5643667
  * [42] Jaroslaw Piwonski and Thomas Slawig  “Metos3D: the Marine Ecosystem Toolkit for Optimization and Simulation in 3-D – Part 1: Simulation Package v0.3.2”  In _Geoscientific Model Development_ 9.10, 2016, pp. 3729–3750  DOI: 10.5194/gmd-9-3729-2016
  * [43] Malte Prieß, Slawomir Koziel and Thomas Slawig  “Surrogate-based optimization of climate model parameters using response correction”  In _Journal of Computational Science_ 2.4, 2011, pp. 335–344  DOI: 10.1016/j.jocs.2011.08.004
  * [44] Malte Prieß, Slawomir Koziel and Thomas Slawig  “Marine Ecosystem Model Calibration through Enhanced Surrogate-Based Optimization”  In _Simulation and Modeling Methodologies, Technologies and Applications_ 197, Advances in Intelligent Systems and Computing  Berlin, Heidelberg: Springer, 2013, pp. 193–208  DOI: 10.1007/978-3-642-34336-0_13
  * [45] Malte Prieß et al.  “Accelerated parameter identification in a 3D marine biogeochemical model using surrogate-based optimization”  In _Ocean Modelling_ 68, 2013, pp. 22–36  DOI: 10.1016/j.ocemod.2013.04.003
  * [46] Nestor V. Queipo et al.  “Surrogate-based analysis and optimization”  In _Progress in Aerospace Sciences_ 41.1, 2005, pp. 1–28  DOI: 10.1016/j.paerosci.2005.02.001
  * [47] Johannes Rückelt et al.  “Parameter optimization and uncertainty analysis in a model of oceanic CO2 uptake using a hybrid algorithm and algorithmic differentiation”  In _Nonlinear Analysis: Real World Applications_ 11.5, 2010, pp. 3993–4009  DOI: 10.1016/j.nonrwa.2010.03.006
  * [48] David E. Rumelhart, Geoffrey E. Hinton and Ronald J. Williams  “Learning representations by back-propagating errors”  In _Nature_ 323.6088 Springer ScienceBusiness Media LLC, 1986, pp. 533–536  DOI: 10.1038/323533a0
  * [49] Jorge Louis Sarmiento and Nicolas Gruber  “Ocean biogeochemical dynamics”  Princeton et al.: Princeton University Press, 2006, pp. 528
  * [50] Charlotte Siberlin and Carl Wunsch  “Oceanic tracer and proxy time scales revisited”  In _Climate of the Past_ 7.1, 2011, pp. 27–39  DOI: 10.5194/cp-7-27-2011
  * [51] T.. Simpson, J.. Poplinski, P.. Koch and J.. Allen  “Metamodels for Computer-based Engineering Design: Survey and recommendations”  In _Engineering with Computers_ 17.2, 2001, pp. 129–150  DOI: 10.1007/PL00007198
  * [52] Alex J. Smola and Bernhard Schölkopf  “A tutorial on support vector regression”  In _Statistics and Computing_ 14.3, 2004, pp. 199–222  DOI: 10.1023/B:STCO.0000035301.49549.88
  * [53] Jacob Søndergaard  “Optimization using surrogate models - by the space mapping technique” Supervisor: Kaj Madsen, 2003, pp. 203  URL: http://www2.compute.dtu.dk/pubdb/pubs/2450-full.html
  * [54] Richard S. Sutton and Andrew G. Barto  “Reinforcement Learning: An Introduction”, Adaptive Computation and Machine Learning series  Cambridge, MA, USA: MIT Press, 2018
  * [55] Pauli Virtanen et al.  “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”  In _Nature Methods_ 17, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2
  * [56] Vito Volterra  “Variations and fluctuations of the number of individuals in animal species living together”  In _Animal Ecology_ New York: McGraw–Hill, 1931, pp. 409–448
  * [57] Carl Wunsch and Patrick Heimbach  “How long to oceanic tracer and proxy equilibrium?”  In _Quaternary Science Reviews_ 27.7, 2008, pp. 637–651  DOI: 10.1016/j.quascirev.2008.01.006
  * [58] Ciyou Zhu, Richard H. Byrd, Peihuang Lu and Jorge Nocedal  “Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale Bound-Constrained Optimization”  In _ACM Trans. Math. Softw._ 23.4 New York, NY, USA: Association for Computing Machinery, 1997, pp. 550–560  DOI: 10.1145/279232.279236