AquaticEcoDynamics / libaed-water

Code for the AED water quality model
GNU General Public License v3.0
6 stars 6 forks source link

Changes to N-fixing penalty function and salinity tolerance #25

Closed MichaelBarryBMT closed 2 years ago

MichaelBarryBMT commented 2 years ago

Hi Matt and Casper

A couple of requests on phyto routines in this one:

Both of these are just suggestions - I just need to know the reasoning behind whatever is correct so that I can finish writing the WQM manual and close this issue out please.

Ta

MB

matthipsey commented 2 years ago

Hi @MichaelBarryBMT

We've done a deep dive into this and made a few changes to clean it up.

In short, the form of the expression is sensitive to S_bep. If S_bep <1 then 0 < f(S) < 1, but if S_bep > 1 then f(S) >1 (except for the salModel 4 which is different).

What Ive done is made a new rule for this :

The user needs to be careful to set S_bep accordingly (either less than or greater than one) or they will get the opposite effect - there is a warning, but this may get missed.

salmodel 4 is a special case, but we've also made a way so that it conforms to the above logic.

I''m not sure how choosing the growth suppression or resp enhancement will affect the calibration process and will need to do some testing to get an idea.

Maybe for your implementation then you could hide one and just rely on a particular form.

https://github.com/AquaticEcoDynamics/libaed-water/commit/712481383d18858414a874cc7ad837fd6b686d7c

aed-modeller commented 2 years ago

A few further notes on the top of Matt's response:

  1. The option "salTol == 1" only works for freshwater, e.g. salinity < S_maxsp;
  2. The option "salTol == 2" or "salTol == 3" work for estuarine water and adapt to a wide range of salinity;
  3. The option "salTol == 4" only works when S_bep < 1

The S_bep is a key parameter for the fSal, when S_bep > 1 then it acts as an amplifier, otherwise is a suppresser.

Below is an example graph showing the response of fSal to S_bep in all salinity limitation functions (note that in future we will use "salTol == -1/-2/-3" for respiration enhancement, and "salTol == 1/2/3/4" for growth suppression):

sal_function

aed-modeller commented 2 years ago

I think we need to add a few lines in the "salTol == 1" option to avoid the fsal goes down below 0. E.g.

IF (salinity>phytos(group)%S_opt .and. salinity<phytos(group)%S_maxsp) THEN
        tmp1 = (phytos(group)%S_bep-1.0) / ((phytos(group)%S_maxsp - phytos(group)%S_opt)**2.0)
        tmp2 = (phytos(group)%S_bep-1.0) * 2.0*phytos(group)%S_opt / &
              ((phytos(group)%S_maxsp-phytos(group)%S_opt)**2.0)
        tmp3 = (phytos(group)%S_bep-1.0) * phytos(group)%S_opt*phytos(group)%S_opt / &
              ((phytos(group)%S_maxsp-phytos(group)%S_opt)**2.0) + 1.0
        fSal = tmp1*(salinity**2.0) - tmp2*salinity + tmp3
     ELSEIF (salinity>=phytos(group)%S_maxsp) THEN
        fSal = 0
     ELSE
        fSal = 1.0
ENDIF
matthipsey commented 2 years ago

Hi @MichaelBarryBMT

Ive had a close look at this.

I think it is implemented correctly in the code as it was, and your suggested modification was not right...

k_nfix is the fractional rate of GPP under full N-fixation (eg 0.8 is 80% efficient growth rate when fixing)

GPP = GPP * ( k_nfix + f(N) [1-k_nfix] )

when f(N) is 1, then the amount in the ( ) sums to 1, and no penalty is applied (as no n-fixing needed)

when f(N) is 0, then the amount in the ( ) sums to k_nfix, and the penalty is fully applied

MichaelBarryBMT commented 2 years ago

Hi Matt Thanks for looking at this. Agreed! MB