Mathematical analysis of the $PO_4$-$DOP$-$Fe$ marine ecosystem model driven by 3-D ocean transport

0.26 MB·English

by  Christina Roschat

Mathematical analysis of the PO -DOP-Fe marine
4
ecosystem model driven by 3-D ocean transport
Christina Roschata,∗, Thomas Slawiga
5
1 aDep. of Computer Science, KielMarine Science - Centre for Interdisciplinary Marine
0 Science, Christian-Albrechts-Universit¨at zu Kiel, 24098 Kiel
2

n
a
J

2 Abstract
2

Marineecosystemmodelsaredevelopedtounderstandandsimulatethebiogeo-
]
P
chemical processes involved in marine ecosystems. Parekh, Follows and Boyle
A
h. introducedthe PO4-DOP-Fe model of the coupled phosphorusand ironcycles
t in2005. Especiallythepartdescribingthephosphoruscycle(PO -DOP model)
a 4
m
is often applied in the context of parameter identification. The mathematical
[
analysispresentedinthis studyisconcernedwiththeexistenceofsolutionsand

1
v thereconstructionofparametersfromgivendata. Bothareimportantquestions
8
2 inthe numericalmodel’s assessmentandvalidationnotansweredsofar. Inthis
4
study, we obtain transient, stationary and periodic solutions (steady annual
5
0
cycles) of the PO -DOP-Fe model equations after a slight change in the equa-
4
.
1
tion modeling iron. This result confirms the validity of the solutions computed
0
5
numerically. Furthermore, we present a calculation showing that four of the
1
: PO -DOP model’s parametersare possibly dependent, i.e. different parameter
v 4
i
X valuesmightbe associatedwith the samemodeloutput. Thereby,we identify a
r relevant source of uncertainty in parameter identification. On the basis of the
a
results,possiblewaysto overcomethis deficitcanbe proposed. Inaddition, the
statedmathematicalconditions forsolvabilityareuniversalandthus applicable
to the analysis of other ecosystem models as well.
Keywords: Marine ecosystem models, Steady annual cycles, Parameter
∗Correspondingauthor. Tel.: [phone]
Email address: [email protected] (ChristinaRoschat)
Preprint submitted toEcological Modelling January 23, 2015
identification
1. Introduction
Being a part of the global carbon cycle marine ecosystems considerably in-
fluence on the earth’s climate. Marine ecosystem models, describing the bio-
geochemical processes involved, provide an important tool to understand these
processes and thus to predict the future concentration of carbon dioxide in
oceans and atmosphere.
One important example for models of this kind is the PO -DOP-Fe model
4
by Parekh et al. (2005), describing the concentrations of three tracers phos-
phate (PO ), dissolved organic phosphorus (DOP) and iron (Fe). Following
4
Kriest et al. (2010), we use the synonymous name N-DOP-Fe model. The
abbreviation N stands for the more general expression “nutrient”.
Fixing the iron concentration, the evolving N-DOP model, consisting of
the first two tracers, describes the marine phosphorus cycle. Despite the low
complexity and associated simplifications, it is still in use (Kriest et al., 2012;
Parekhet al., 2006). An assessment on the basis of oceanic observations even
indicates that, under certain circumstances, the N-DOP model can compete
with more complicated ones (Kriest et al., 2010, 2012). Due to the low com-
plexity, it is often used for the purpose of testing numerical methods, e.g. in
Prieß et al. (2013).
Thethirdtracerexpressestheinfluenceofirononthephosphoruscycle. Iron
is an important nutrient. For instance, Kriest et al. (2012) assume that their
observedmisfitofmodeloutputandobservationaldatamightbeduetomissing
ironlimitation. Themodelcanalsopicturethechangesinthephosphoruscycle
induced by artificial fertilization with iron.
The N-DOP-Fe model consists of three advection-diffusion-reaction equa-
tions,eachcharacterizingonetracerconcentrationonathree-dimensionalocean
domain. The concentrations are influenced by ocean transport, i.e. advection
anddiffusion, andbiogeochemicalprocessesmodeledbyspecific reactionterms.
2
Thefirstpartofthisworkstudiestheexistenceoftransient,periodicandsta-
tionarysolutionsoftheN-DOP-Femodelequations. Periodicsolutions(steady
annualcycles)arethemostrelevantones. Suchresultsaredesirablebecausethe
numerical model output used in applications is computed by approximating a
solutionoftheoriginalmodelequations. Therefore,thequalityandreliabilityof
the model output depends on the solvability of these equations. A comparable
analysis of the N-DOP-Fe model equations has not been undertaken so far.
Any model has to be calibrated, i.e. adapted to the ecosystem in ques-
tion. This is achieved by the choice of the model’s parameters. Parameters,
like e.g. remineralization rates or half saturation constants, are essential quan-
tities characterizing the processes modeled. The N-DOP model includes seven
parameters.
Parameters can be estimated by means of laboratory experiments. If such
measurementsaretoodifficultorexpensive,parametersarealternativelyidenti-
fied via optimization with respect to observational data. Here, a parameter set
minimizing the difference between model output and data is determined. An
important area of application of the N-DOP model is the testing of numeri-
calmethods designedto solvesuchkindofminimizationproblems(Prieß et al.,
2013). The method in question is applied to synthetic data, i.e. data corre-
spondingto knownoptimalparameters,assumingthatacorrectmethodisable
to identify these. However, this assumption only holds true if all parameters
are uniquely identifiable, i.e. if each possible model output is associated with a
single parameter set.
The second part of this work is dedicated to the question, unanswered so
far,whichoftheN-DOP modelparametersareuniquelyidentifiable. Wename
and justify different alternatives to alter the reaction terms in such a way that
all parameters become identifiable.
The mathematical results about the N-DOP-Fe model are introduced and
explainedto readersfromother disciplines than mathematics. By outlining the
mathematicalproceeding,we intend to provideanoverviewaboutwhichmath-
ematicalconditions andassumptionsareresponsiblefor the model’s properties.
3
Thesemaybeaguidelinefortheinvestigationanddevelopmentofothermodels
or alternative reaction terms.
The paper is structured as follows: In the next section, we introduce the
classicalandthe weak formulationofthe N-DOP-Femodel equations,the lat-
terbeingtheobjectofinvestigation. Thenextthreesectionseachdealwithone
type of solution. The corresponding results concerning existence and, where
possible, uniqueness are formulated and shortly justified. In Sec. 6, we inves-
tigate identifiability of parameters. In the last two sections, the results are
discussed and conclusions are drawn.
2. Model equations
ThefollowingintroductionoftheN-DOP-Femodelisbasedontheoriginal
work by Parekh et al. (2005) as well as on Roschat and Slawig (2014a,b).
The ecosystem is located in a three-dimensional bounded domain Ω deter-
mined by the bounded water surface Ω′ ⊆R2 and the depth h(x′)>0 at every
′ ′ ′
surface point x ∈ Ω. The boundary Γ is the union of the surface Γ and the
boundary inside the water.
Thedomainisseparatedintotwolayers,the euphotic,light-floodedzoneΩ
1
uptoadepthofh¯ :=120mandthedark,aphoticzoneΩ beneath. Theactual
e 2
′ ′
depthoftheeuphoticzonebeneathsomesurfacepointx isdefinedbyh (x):=
e
min{h¯ ,h(x′)}. We split the surface into the part Ω′ := {x′ ∈ Ω′;h(x′) > h¯ }
e 2 e
′ ′ ′
abovetheaphoticzoneandtherestΩ :=Ω \Ω . Theboundaryisanalogously
1 2
divided into the euphotic part Γ and the aphotic part Γ .
1 2
For the sake of a clearer terminology, we write henceforth y for the phos-
1
phate concentration, y for the concentration of DOP and y for the iron con-
2 3
centration. The three tracers, assembled in the vector y = (y ,y ,y ), are
1 2 3
characterizedby the system of advection-diffusion-reactionequations
∂ y +div(vy )−div(κ∇y )−λy +d (y) = 0
t 1 1 1 2 1
∂ y +div(vy )−div(κ∇y )+λy +d (y) = 0 (1)
t 2 2 2 2 2
∂ y +div(vy )−div(κ∇y )+J −λy R +d (y) = S ,
t 3 3 3 Fe 2 Fe 3 Fe
4
eachdependentonspaceandtimeinafinitetimespan[0,T]. Regardingperiodic
solutions, a reasonable value for the final point of time T is one year.
Thesecondandthirdtermsontheleft-handsidesdescribetheoceandynam-
ics by means of the velocity field v and the diffusion coefficient κ. Sometimes
these terms are summarized in a linear operator or, if already discretized, in a
matrix (Kriest et al., 2010, 2012). To reduce computational effort, the models
are often run in an “offline” mode, i.e. the influence the tracers have on the
oceans dynamics is neglected. This has the advantage that the values of v and
κ can be precomputed with one run of an ocean circulation model. For that
reason, v and κ are assumed to be known in the theoretical analysis as well.
Since turbulent exceeds molecular diffusion by far the values of κ are assumed
equal in all equations.
Furthermore, since the ocean is a closed system, it is reasonable to assume
that the velocity field v is divergencefree anddoes notpoint out of the bound-
aries. These two properties are crucial for the mathematical analysis.
The reactionterms d (y),J ,λy and λy R describe the biogeochemical
j Fe 2 2 Fe
coupling. Therefore, they depend on one or several of the tracers. On the
right-hand side of each equation, sources and sinks of the respective tracer are
displayed. Only iron has a non-zero source term.
Inthe followingtwosubsections,weintroduceandderivethe reactionterms
associated with the phosphorus cycle and the iron equation, respectively. Two
further subsections deal with boundary conditions and weak solutions.
2.1. The phosphorus cycle
One important process is the remineralizationof y into y . In the first two
2 1
equationsofEq.(1),it ismodeledas afirstorderlossprocesswitha remineral-
izationrate λ between 0 and 1. Being independent of light, this transformation
takes place in the whole domain Ω. The remaining processes, represented by
the reaction terms d and d , differ according to the layers. In the euphotic
1 2
zone, y is consumed by the photosynthesis of phytoplankton. This uptake is
1
5
modeled by
y y Ie−x3KW
1 3
G(y ,y ):=α .
1 3 |y1|+KP |y3|+KF |Ie−x3KW|+KI
Here, the maximum uptake α > 0 is limited by the present concentrations of
phosphate y and iron y as well as insolation via Michaelis-Menten kinetics.
1 3
The corresponding saturation functions are equipped with the half saturation
constantsK ,K ,K . Theabsolutevaluesinthedenominatorsensurethatthe
P F I
fractions are mathematically well-defined because it is a priori unknown if the
solutions y and y are nonnegative. If this holds true, the above formulation
1 3
is in accordance with Michaelis-Menten kinetics. Incidence of light is formal-
′
ized by the insolation I = I(x,t) depending on the water surface and time.
The incident light decreases exponentially with depth x and the attenuation
3
coefficient K for seawater. Below the euphotic zone, the values of I are zero
W
(Paltridge and Platt, 1976).
A fraction ν ∈ (0,1] of the uptake G is transformed into y while the rem-
2
nants, integrated over the whole water column, are exported into the aphotic
zone Ω . The remineralization during the sinking of particles is described by a
2
parameter b.
The outlined processes are represented by nonlinear coupling terms with
different appearance in each layer. In the euphotic zone Ω , they are given by
1
d (y):=G(y ,y )
1 1 3
d (y):=−νG(y ,y ),
2 1 3
and, in Ω , by
2
he b x −b−1
d (y):=−(1−ν) G(y ,y )dx 3
1 1 3 3h¯ ¯h
Z0 e (cid:18) e(cid:19)
d (y):=0.
2
We remark that the N-DOP model contains the seven parameters λ,α,K ,
P
K ,K ,b,ν.
I W
6
2.2. The iron equation
In opposite to phosphorus, there is a source term for iron given by S :=
Fe
βF which is non-zero only in the euphotic zone. The parameter β represents
in
the solubility of iron in seawater and F quantifies the aeolian source of iron.
in
The reaction terms express how the phosphorus cycle, complexation and
scavenginginfluenceontheironcycle. Theironconcentrationincreaseswithre-
mineralizationanddecreaseswithconsumptionofphosphorus. Beingexpressed
in phosphorus units, these values are multiplied by the constant ratio R in
Fe
order to become iron units. Thus, iron increases with λy R . The decrease
2 Fe
induced by consumption is expressed by the coupling term
d (y)=G(y ,y )R .
3 1 3 Fe
2.2.1. Scavenging and complexation - original formulation
The summand J represents the influence of complexation and scavenging
Fe
on the iron concentration. Since a certain amount of iron is complexed with
′
organic ligand, the total iron y is split into free iron Fe and complexed iron
3
′
FeL. Similarly, the total ligand L is the sum of free ligand L and the com-
T
plexed ligand, equal to the complexed iron FeL. These relations are expressed
′ ′
by the formulae y = Fe +FeL and L =L +FeL. Parekhet al. (2005) set
3 T
L =1.
T
Onlythe freeironissubjectto scavengingwhichisthusmodeledasthefirst
order loss process J := τk CΦFe′. Here, τ,Φ are numbers, k is the initial
Fe 0 p 0
scavenging rate and C represents the particle concentration which decreases
p
with depth.
Since a feasible mathematical formulation requires that J depends at
Fe
′
least on one of the tracers y ,y or y , we express Fe using y . To this
1 2 3 3
end,Parekh et al.(2005)additionallyprovidethe equilibriumrelationshipK =
′ ′
FeL/FeL with a positive constant K. Inserting the equivalent expression
′ ′
L =FeL/KFe into the equation for ligand we obtain
FeL 1
′
L =FeL+L =FeL+ =FeL 1+ .
T KFe′ KFe′
(cid:18) (cid:19)
7
′
FeL can be replaced by y −Fe. This gives
3
1 y 1
′ ′ 3
L =(y −Fe) 1+ =y −Fe + − .
T 3 KFe′ 3 KFe′ K
(cid:18) (cid:19)
With the abbreviation H(y ):=L +1/K−y , this proves equivalent to
3 T 3
y
Fe′2+H(y )Fe′− 3 =0.
3
K
The two solutions of this quadratic equationare known. It can be easily shown
that one of them has only negative values and is therefore inappropriate to
describe the amount of free iron. Thus, we find
1 H(y )2 y
′ 3 3
Fe(y )=− H(y )+ + .
3 3
2 4 K
r
To ensure that the square root is real we show that the radicand r := (L +
T
1/K−y )2/4+y /K is positive. This is obvious whenever y is positive, since
3 3 3
in this case r ≥y /K >0. In the non-positive case, we transform r into
3
2
1 1 y
3
r= L + −y +
T 3
4 K K
(cid:18) (cid:19)
2
1 1 1 1 1 y
= L + − L + y + y2+ 3
4 T K 2 T K 3 4 3 K
(cid:18) (cid:19) (cid:18) (cid:19)
1 1 2 1 1 1
= L + + y2− y L − .
4 T K 4 3 2 3 T K
(cid:18) (cid:19) (cid:18) (cid:19)
SinceKtendstobealargenumberandL takesavaluearoundone(Parekh et al.
T
(2005) set K = exp(11) and L = 1) the expression L −1/K is nonnega-
T T
tive. Therefore, y ≤ 0 yields −y (L −1/K)/2 ≥ 0 and, as a consequence,
3 3 T
r ≥(L +1/K)2/4>0. In particular, this result implicates that Fe′ is differ-
T
entiable everywhere on R. It can be shown that the derivative is positive and
′
thus that Fe is a monotonically increasing real function.
In Fig. 1, we see that, as long as approximately y <L =1, the values for
3 T
′
Fe increase very slowly, for y > L they increase with a gradient of almost
1 T
one. This behavior corresponds to the statement of Parekhet al. (2005) that
′
they “rapidlyprecipitate Fe when Fe >L ” meaning that, as long as ligand
T T
is available, only a small amount of iron remains free.
8
0.6
0.5
0.4
0.3
Free iron 0.2
0.1
0
−0.1
−0.5 0 0.5 1 1.5
Total available iron
Figure1: GraphoffreeironFe′independenceoftotalavailableirony3assumingK=exp(11)
andLT =1(Parekhetal.,2005).
Thus, the reaction term associated with scavenging and complexation de-
pending on the third tracer y is determined by
3
J (y ):=τk CΦFe′(y ).
Fe 3 0 p 3
2.2.2. Scavenging and complexation - adjusted formulation
The reactiontermJ (y ) derivedabovelacksone importantpropertywith
Fe 3
regardto solvability. This is due to the fact that the function
1 1 1 1 y
′ 3
Fe(y )=− L + −y + L + −y +
3 T 3 T 3
2 K s4 K K
(cid:18) (cid:19) (cid:18) (cid:19)
modeling free iron levels off fast in negative direction. For that reason, we
′
propose an alternative term FeF(y ) for free iron. Fig. 1 indicates that Fe
3
resemblesastraightlinewithaslopeof1forapproximatelyy >L . Ify ≤L ,
3 T 3 T
′
the curve tends to zero very fast. As a substitute for Fe we therefore define a
piecewise linear function, composed of a line with a slope of 1 and another line
′
through zero. As their “meeting point” we determine (L ,Fe(L )) where the
T T
′
function Fe apparently turns upwards. Formally, the function FeF depending
9
0.03
0.025
0.02
Free iron 0.015
0.01
0.005
0
−0.005
0.94 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04
Total available iron
Figure 2: Graph of free iron Fe′ (red line) in comparison to the alternative function FeF
(blueline)aroundy3=LT assumingK=exp(11)andLT =1(Parekhetal.,2005).
on total iron y is defined by
3
′
y +Fe(L )−L if y >L ,
3 T T 3 T
FeF(y )=
3 
FeL′(TLT)y3 if y3 ≤LT.

′
Figure 2 shows the difference between exemplary curves of Fe and FeF
around y = L . Comparing the complete curves reveals that, in the positive
3 T
region,thepiecewiselinearfunctionFeF liesslightlyabovetheoriginalandtheir
distance remains small. In the negative region, the curves behave conversely.
Thedistanceincreaseswithdecreasingvaluesofy becauseFeF hasaconstant
3
′
slopeandFe adecreasingone. However,negativevaluesfory arenotrelevant
3
in most applications.
Takingintoaccounttheseconsiderations,apossiblealternativeforJ might
Fe
be the adjusted reaction term
J(y ):=τk CΦFeF(y ).
3 0 p 3
10