Fix GR pathology

Issue initially reported to Vazken and Charles by Wouter Knoben on the Excel GR4J in 2017, also reported recently by @laurent.strohmenger on GR5J recently.

Note: this issue concerns airGR but is not specific to the package, it comes from the GR implementation.

Under some conditions (i.e. some parameters sets, some climatic conditions), the lower part of the GR models, namely the routing store and the exchange component, enter in a mode where it cannot escape anymore. Due to a positive X2 and a positive exchange larger than the routing store outflow, the routing store comes to an equilibrium and a non null discharge is produced, even though no more rainfall is injected.

This is very annoying, as the simulation quality becomes very poor. This issue may not be found under calibration, because a low objective function would lead to avoid these parameter sets, but some parameter sets can look ok under calibration and then give erroneous simulations, under a different climate for instance (e.g. high rainfall).

Included is the email exchange with Wouter Knoben colleagues.

GR4J_pathology_Knoben.docx

Attached also an R code reproducing the bug proposed by @laurent.strohmenger and augmented by the translation in R of the GR5J routing part producing the bug. It shows that X2 does not need to be larger than X3, oppositely to what was discussed in 2017.

GR5J_bug.R

There is no solution to date, but this must be cured.

added BUG FORTRAN labels

Here is a function that identifies potentially toxic parameter sets. The idea behind the math is that Exch should never be greater than QR. Thus, Exch=f(Rout) should be less than QR=f(Rout) for Rout where dExch/dRout == dQR/dRout.

check_params <- function(x2, x3, x5){ 
  #return TRUE for toxic parameters set
  dExch <- x2/x3 
  S <- ((dExch * (4*x3^4))/5)^(1/4) #Rout when dExch/dRout == dQR/dRout
  
  Ex_tan <- x2 * (S/x3 - x5) #Exch for Rout = S
  QR_tan <- S*(1.-(1.+(S/x3)**4.)**(-1./4.)) #QR for Rout = S 
  
  if (is.na(S) | Ex_tan < QR_tan){
    return(FALSE)
  }
  return(TRUE)
}

*Edit: The derivative dQR/dRout is approximative, this is not a perfect solution.

@laurent.strohmenger

In the GR5J code, I've made a small modification: I check the difference between the exchange and the outflow from the routing store. If this difference becomes higher than -0.1 for 5 consecutive time steps (I tried with strictly positive, but the asymptote can be going upwards, so it did not work), I remove 1 mm in the routing store. This seems:

• to allow preventing the impact for toxic parameters
• not to change a working well simulation

Could you please extend the tests since you have a large set of working and non-working catchments? You'll have to rebuild an airGR package to test it of course.

frun_GR5J.f90

@guillaume.thirel The modification seems to allow the model to escape to the loop for the 2 stations that were affected by this problem.

After some quick tests, calibrating the "healthy" catchment with the modified GR5J caused substantial changes in the parameter values.

changed title from GR pathology to Fix GR pathology

added ENHANCEMENT label

changed milestone to %v2.0