a possible bug in fennel biological model

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
tony1230
Posts: 87
Joined: Wed Mar 31, 2010 3:29 pm
Location: SKLEC,ECNU,Shanghai,China

a possible bug in fennel biological model

#1 Unread post by tony1230 »

I cannot confirm absolutely that this is a bug report, but i find something...

According to Katja's paper/2006, the grazing term on Phytoplankton by Zoop is defined as g*Zoo(page 3), where g is the grazing rate of Phyt. And g is represented by Holling-type s-shaped curve as(refer to P.4)

Code: Select all

g = gmax*Phyt*Phyt/(kp + Phyt*Phyt)
,thus, the portion of Phyt grazed by Zoop(i call 'Zphy' here) should be

Code: Select all

Zphy = g*Zoo = gmax*Zoo*Phyt*Phyt/(kp + Phyt*Phyt)
But when i looking into the file fennel.h (Nonlinear/Biology/),i found this

Code: Select all

! Phytoplankton grazing by zooplankton.
!
              cff1=fac1*Bio(i,k,iZoop)*Bio(i,k,iPhyt)/                  &
     &             (K_Phy(ng)+Bio(i,k,iPhyt)*Bio(i,k,iPhyt))
              cff3=1.0_r8/(1.0_r8+cff1)
              Bio(i,k,iPhyt)=cff3*Bio(i,k,iPhyt)
.
As is calculated from here, Zphy = Bio(i,k,iPhyt) = cff3 * Bio(i,k,iPhyt) not equal to above formula. The term cff3 here should be the cff1,i.e. Bio(i,k,iPhyt)=cff1*Bio(i,k,iPhyt). IF this is the case, can that make the sense.

- Shou

often

Re: a possible bug in fennel biological model

#2 Unread post by often »

In fennel.h, implicit algorithm is used to prevent negative value.

explicit algorithm

Code: Select all

Zphy = g*Zoo = gmax*Zoo*Phyt*Phyt/(kp + Phyt*Phyt)
Phyt(new)=Phyt(old)-Zphy
Zoo(new)=Zoo(old)+Zphy
implicit algorithm

Code: Select all

cff1=gmax*Zoo*Phyt/(kp + Phyt*Phyt)
cff3=1/(1+cff1)
Phyt(new)=Phyt(old)*cff3
Zoo(new)=Zoo(old)+Phyt(new)*cff1
cff1 look like missing an *Phyt, but I don't think that's a mistake.

User avatar
shchepet
Posts: 188
Joined: Fri Nov 14, 2003 4:57 pm

Re: a possible bug in fennel biological model

#3 Unread post by shchepet »

There is no bug here as far as I can see.

Code: Select all

cff1=gmax*Zoo*Phyt/(kp + Phyt*Phyt)
cff3=1/(1+cff1)
Phyt(new)=Phyt(old)*cff3
Zoo(new)=Zoo(old)+Phyt(new)*cff1
The whole biological time stepping procedure is dated back into 1998 -- about 14 years ago, see http://www.atmos.ucla.edu/~alex/biology, and the main idea is splitting by physical processes with implicit treatment of the component being consumed.

It is sometimes known as method of fractional steps.

An obvious motivation to use implicit step, is that the resultant Phyt(new) would never be negative. However, there are other reasons as well. Three things needs to be explained to understand why it is chosen to be this way:
  • (i) cff1 is proportional to time step "dt" (it is not self-obvious from the code fragments quoted above, but people should be able to trace back), so if time step is sufficiently small, so does cff1, so cff3=1/(1+cff1) can be expanded in Taylor series as 1-cff1, which corresponds to the explicit forward Euler stepping. So the two versions, explicit and implicit, converge to the same limit when dt --> 0.

    (ii) It can be immediately verified that Phyt(new)+Zoo(new)=Phyt(old)+Zoo(old), no matter what, so the procedure is conservative for the net content. (This is why the last line of the code, Zoo(new)=Zoo(old)+Phyt(new)*cff1 has Phyt(new) in the r.h.s, rather than the old-time-step value.

    (iii) The mathematical nature of the system (biological equations) is not just relaxation toward equilibrium (as would be in a system of parabolic type), but may possess more complex (actually quite rich) behavior: oscillations about phase equilibrium, limits cycles, loss of stability of a stationary point, etc... As forward Euler step is unconditionally unstable for an oscillatory system, and biological system may possess such behavior, ignoring this leads to physically wrong solutions (some major biological codes are ignoring this basic fact). The procedure consisting of partially implicit fractional steps is deliberately designed biased toward non-oscillation (and actually numerically induced dissipation), so any oscillations (usually seen as flashes of concentrations of components) are most likely physical, and not artifacts of numerical instability.

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Re: a possible bug in fennel biological model

#4 Unread post by arango »

No, there is not a bug here. This is part of the implicit algorithm. Sasha explained this well above. All the biological models in ROMS have extensive comments about this. You just need to follow the comments and algebra between producers and consumers.

This is not the first time that we get this type of messages about the ecosystem models. The logic of these algorithms is not that difficult to follow. Everything is explained well.

Post Reply