a possible bug in fennel biological model

[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]

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.

[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.

[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.