Fasham's Biological Model as a Dynamical System:
Overview of Forced Regimes
The biological module of ROMS was updated in order to make it more robust.
Now the time stepping is done implicitly (backward Euler step).
The design goal here is robustness, as well as non-oscillatory behaviour
of the numerical solution.
In the case when natural decay rates (mortality and sinking of Zoo and
Phyto, as well as loss of Detritus to sinking) are set to zero, the
numerical procedure conservs the sum of the all five components exactly
on the time-discretized level.
Since the system of differential equation is nonlinear, and iterative
approach is used to achieve the implicit solution.
During each iteration a sequence of sub-steps is performed to update
the concentration of each biological component in a such a way that the
just computed values are immediately used in computation of the next
component(approach sometimes known as splitting by physical
precesses ).
The sequence of sub-steps is according to the main food chain,
NO3 --> NH4 --> Phyto --> Zoo --> Det --> next iteration .
The component being consumed is treated implicitly during the sub-step
in which it is consiumed, which guarantees positive definiteness,
no matter how fast the consumption rate is or how large is the time step.
The new file bilogy.F is available here .
It is compatible as is with the code in Scripps.
Forced Seasonal Cycle
In the following series of experiments the biological model was extracted
from ROMS and was run as a box model, i.e dynamical system.
Technically it is done by running the same code, where all model
arrays are declared (1,1,1).
The model was forced by applying seasonally varying source of NO3:
source_NO3 = A * (1.-cos(2*pi*tdays/360)
where A is supply of NO3 dimensioned as [mMol N /(m^3 day)].
In addition to that the model is parametrically forced through solar radiation,
which has dayly cycle. This forcing does not affect directly any biological
component (in sense of increase or decrease of concentration), but regulates
the transformation rate of NO3 and NH4 into Phytoplankton.
In the case of no external NO3 supply the system goes to the state of zero
concentrations of Phyto- and Zoo-plankton and Detritus (according to the
e-folding scales of their decay), with NO3 consumed first, and some
nonzero concentration of NH4 remains after everything else disappears.
In the case of periodic forcing the system reaches seasonal cycle after
only a couple of years for moderate values of A. It takes longer,
when the supply is small because the initial condition (all concentrations
are set to 0.1 mMol N/m^3 is too far from equilibrium.
Fig. 1 A typical multi-year history for the transitional
(in sense defined later) forcing amplitude A is where the cycle repeats
itself almost exactly, including small details already after 2 years.
Amplitude of NO3 source A = 0.0425. First two years of this
experiment will be shown in more details on Fig. 10 below.
Fig. 2 For a very small NO3 flux it takes several years to reach
a periodic regime and the adaptation is not straightforward
(20 year history is shown here, A = 0.00002). It should be noted
that in this setup there is always surviving Phytoplankton, no matter how small
the supply of NO3 is: since NO3 is stable, and since the NO3 uptake is
proportional to the product of concentrations of NO3 and Phyto (in the limit
when Phyto concentration is small), the system accumulates enough NO3 to make
the uptake balance the mortality and sinking of Phyto.
In contrast to that, survival of Zooplankton is possible only when the
concentration of Phyto exceeds certain threshold, which is possible only
when there is a sufficient supply of NO3. When no Zooplankton is present,
within wide range of parameters, in time-averaged sense, the governing
balance is,
[Phyto] * [NO3] ~ source_NO3
and, moreover, increase of source of NO3 does not increase concentration of
NO3, but rather goes to Phytoplankton. In the two cases shown above the
equilibrium concentration of NO3 is around 0.03 mMol N /m^3, despite the fact
that supply of NO3 differs by three orders of magnitude.
For most of the parameters studied below 2 years is representative enough
(longer runs were also made to verify this statement, but not shown here).
First transition: Survival
of Zooplankton
Fig. 3 The source amplitude is 0.001. The duration of this run
is 2 years. The source of NO3 has maximum in the middle of the year, and zero
at the beginning - end (vertical dashed lines).
Fig. 4 The source amplitude is 0.01.
This image and image above are typical
seasonal cycles for small supply. Note the relative positions of Green
and Blue curves. Despite the supply of NO3 was increased by one order of
magnitude, the equilibrium concentrations of NO3 and NH4 here are within
0.01 ... 0.03 range. Blue and violet curves look like wide bands.
This is due to dayly cycle.
In both cases Zooplankton irreversibly disappears.
Fig. 5 The source amplitude is 0.015.
Zooplankton marginally survives,
or close to survival, since the concentration
of Phytoplankton exceeds the critical value at some interval of time.
In this regime Zooplankton concentration is too small to control
anything. Note that the shape of the green curve is virtually the same as
in the case above.
Fig. 6
The source amplitude is 0.02, that is only 1.5 times larger than in the
preceding case. Peak Zooplankton concentration has increased by more than
order of magnitude and now the concentration of Phyto is contained.
Note deformation of the shape of green curve.
Fig. 7
The source amplitude is 0.03.
Second transition:
Appearance of 30-day oscillations.
Fig. 8
The source amplitude is 0.035.
Zooplankton is strong and controlling.
No 30-day (or similar period) oscillations are present thus far.
Fig. 9
The source amplitude is 0.04. Small oscillations are present.
Fig. 10
The source amplitude is 0.0425. This regime is
transitional to oscillatory regime. Most oscillating components here are
Phytoplankton and NH4. Zoopkankton (with sufficient concentration)
mostly plays catalytic role here.
One should note that there are two Phyto - NH4 cycles in this system:
NH4 --> Phyto(mortality) --> Det(decay) --> NH4 (slow)
NH4 --> Phyto(grazing)--> Zoo(excretion)--> NH4 (fast)
With the rate of the second one proportional to Zoo concentration.
However, since most of Zoo grazing (75%) results into conversion of Phyto into
Det and NH3, Zoo concentration itself does not oscillate very much.
Fig. 11
The source amplitude is 0.045
Fig. 12
The source amplitude is 0.05
Fig. 13
The source amplitude is 0.06
Note that the sum of all five components, the black curve above, does not
oscillate.
Third transition: Oscillation
dominated (non-physical?) regime
Fig. 14 The source amplitude is 0.08.
Fig. 15 The source amplitude is 0.1.
Fig. 16 The source amplitude is 0.2
Fig. 17 The source amplitude is 0.4
Fig. 18 The source amplitude is 0.8.
Looks ugly, yeah?. What one can see here:
Phyto grows, until it exceeds critical value for grows of Zooplankton.
Once Zoo grows, it happened so fast that Phyto is consumed instantly.
Nothing to graze, Zoo decays at its natural mortality-sinking rate,
so does Det (at a slower rate). Unconsumed NO3 and NH4 go to very high
values, providing environment for explosive growth of Phyto in the
next cycle.
Regime transition of the system
without dayly cycle of incoming solar radiation
In this series of experiments we have replaced the dayly cycle
modulation of PAR with its equivalent average value (which is
1/pi). Overall, this change shifts the threshold of oscillations
toward somewhat larger values of NO3 source (A=0.055 vs.
0.0425), but overall behaviour is similar to what we saw above.
Here is transition toward the oscillatory regime (diagrams with
small values of source amplitude are similar to that above, and
therefore not shown).
Fig. 19 The source amplitude is 0.045.
There is no evidence of oscillations thus far.
Fig. 20 The source amplitude is 0.05.
Small oscillations are present.
Fig. 21 The source amplitude is 0.055.
Fig. 22 The source amplitude is 0.06.
Fig. 23 The source amplitude is 0.1.
Fig. 24 The source amplitude is 0.2.
Regime transition of the
modified system (1)--(5) as suggested by John Moisan
As suggested by John, the denomenator in the Zooplankton term,
see Eq.(8) above, was removed:
G = g_max * Phyto * Zoo
instead of the original
G = g_max * Phyto * Zoo / ( K_Phyto + Phyto)
(it should be noted that in the Fasham's choice of parameters
K_Phyto= 1 mMol N /m^3, so that its absence in denominator is
OK, as long as the concentration is measured in [ mMol N /m^3]).
This measure has two concequences:
(1) it increases overall grazing of Phytoplankton (thought
insignificantly for very small concentrations of Phyto),
and,
(2) it removes the saturation effect on the Zooplankton grazing
in the case of large concentration of Phyto.
Thought it was claimed that it removes the 30-day oscillations,
I found that it just shifts the threshold of the oscillatory regime
toward larger NO3 suppy.
Fig. 25 The source amplitude is A=0.04.
This should be directly comparable to the case shown on
Fig. 9, which has the same
NO3 supply and differs only by the absense of saturation in the
grazing term.
The obvious difference from the case on Fig. 9 is the significant
suppression of the spring bloom of Phytoplankton.
Case shown on Fig. 9 was transitional toward 30-day oscillations:
minor increase in supply triggres them.
No oscillations are seen in the present case, neither in the next
one.
Fig. 26 The source amplitude is A=0.05.
Note that the Phytoplankton curve is almost flat: there is no
spring bloom in the second year.
No ocsillations thus far.
Fig. 27 The source amplitude is A=0.06.
Transition to oscillatory regime, thought the amplitude of
oscillations is somewhat smaller than that on
Fig. 13, which has
the same supply of NO3.
Fig. 28 The source amplitude is A=0.07.
30-day oscillations are developed.
Response to constant
and slowly varying flux of NO3: Stationary points and loss
of stability
In the case when supply NO3 is constant in time the system
Eq (1)-(5) has stationary solutions, which may be both stable
and unstable, depending on the level of supply.
Fig. 29
Diagram of stationary fixed points of system Eqs. (1)-(5).
Parameters changing with dayly and seasonal cycles are
replaced with their time-averaged values.
Horizontal axis: NO3 source from 0.0 to 0.008,
linear scale.
If NO3 supply is below approximately 0.0015, Zooplankton
cannot survive. In this regime Phytoplankton concentration
is determined by the balance of NO3 supply and mortality of
Phytoplankton. Note that concentration of NO3 practically
does not depend on its supply within this range.
Once the supply is above 0.0015, Phytoplankton concentration
is controlled by Zooplankton and remains constant with the
increase of supply.
The fact that Phyto is constant in this case is obvious
from the structure of r.h.s of Zooplankton Eq. (5).
Both the decay term and the grazing term there are linearly
proportional to the concentration of Zoo.
Thus, the r.h.s. of Eq. (5) vanishes if
[Zoo] = 0
or if
mu_[Zoo] +
mu_[Zoo][NH4] +
mu_[Zoo][Det] =
beta g_max [Phyt] /
(K_Phyt+[Phyt])
The first one, formally speaking, remains stationary point
regardless of the value of NO3 supply, however, when supply
exceeds the first critical value of 0.0015, it becomes
unstable.
Above this critical value the second equation determines
the stationary concentration of Phytoplankton.
Once the supply reaches critical value of 0.0077, the
stationary solution becames unstable.
Fig. 30 Temporal response of the system to the
quasi-stationary supply.
Here NO3 supply linearly changes from 0.0075
to 0.0079 over the period of 40 years.
After the initial adjustment, which takes approximately 1/2
years, the sustem reaches equilibrium state, which is slowly
changes according to the change of supply.
Once it passes the critical value of approximately 0.0077,
the stationary state becomes unstable.
Fig. 31 Similar to Fig. 30, but the supply is
varying from hi to low: from 0.0079 to 0.0075.
If compared to Fig. 30, note the effect of hysteresis:
transition from non-oscillatory to oscillatory regime
happens at approximately 0.0077, while backward transition
at 0.0075.
Fig. 32 Phase trajectories corresponding to the
limit cycle of the system with constant supply.
This integration was started from the instable stationary point.
NO3 supply here is 0.00765, which is
approximately 0.1% greater than the critical value of 0.007643,
corresponding to the loss of stability.
Note almost perdect elliptical shape of each curve.
Fig. 33 Phase diagram similar to that on Fig 32,
but for slightly larger NO3 supply 0.008.
Fig. 34 Similar to Fig 32, but NO3 supply 0.016
Fig. 35 Similar to Fig 32, but NO3 supply 0.025