Opened 4 years ago

Closed 4 years ago

#412 closed defect (fixed)

Very likely bug in the calculation of k (hydraulic conductivity)

Reported by: aducharne Owned by: aducharne
Priority: critical Milestone: ORCHIDEE 2.0
Component: Physical processes Version: trunc
Keywords: Cc:

Description

La conductivité hydraulique k varie dans le modèle en fonction de ks et de mc, qui varient tous deux en fonction de la profondeur donc de la couche jsl.

La valeur de k qui est utilisée dans hydrol_soil, notée k(ji,jsl), est définie dans hydrol_soil_coef pour chaque soiltile jst, à partir des valeurs a_lin et b_lin qui permettent de linéariser K(theta) en 50 segments pour résoudre l'équation de Richards, via une matrice tridiagonale.

hydrol_soil_coef sort en fait 4 variables a, b, d, k, qui sont essentielles pour l'hydrodynamique, et non indépendantes:

  • k=a*mc+b (en principe) et d=f(k)
  • a et d servent dans hydrol_soil_setup pour définir des paramètres ensuite intégrés dans les matrices traitées par tridiag
  • b sert à construire les matrices traitées par tridiag
  • k sert dans hydrol_soil_infilt (calcul de l'infiltration et première actualisation de mc), puis pour le calcul du drainage en cohérence avec le système tridiag.

Une autre chose que fait hydrol_soil_coef est de modifier k et les variables associées par kfact_root qui augmente k vers al surface du sol.

Mais il semble y avoir une incohérence entre les valeurs de k d'une part et a,b,d d'autre part qui sortent de hydrol_soil_coef, à cause du facteur kfact_root. C'est vrai dans les deux cas ok_freeze_cwrr true ou false, même si je donne ici que le cas sans gel:

DO jsl=1,nslm
          DO ji=1,kjpindex 
             
             ! it is impossible to consider a mc<mcr for the binning
             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji)))
             
             i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin)
             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
             k(ji,jsl) = MAX(k_lin(imin+1,jsl,njsc(ji)), &
                  a_lin(i,jsl,njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji)))  ! in mm/d
          END DO 
       END DO

Comme a_lin, b_lin et k_lin (calculés au préalable dans hydrol_var_init ne dépendant pas de kfact_root, on a donc, il me semble, une sous-estimation de k du facteur kfact_root (lequel sert à augmenter k).

Ceci doit donc entraîner une sous-estimation de l'infiltration et du drainage, ce dernier d'un facteur plus faible car kfact_root diminue avec la profondeur. Quelles sont les conséquences sur le bilan d'eau, l'ET et la croissance des plantes ??? Doit-on corriger ce bug pour CMIP6 ?

La solution est simple, mais je pars en réunion... je compléterai plus tard.

Attachments (3)

bug_ksat_Nov2017.pdf (68.2 KB) - added by aducharne 4 years ago.
hydrol_bugksat.f90 (361.7 KB) - added by aducharne 4 years ago.
hydrol_bugksat2.f90 (362.1 KB) - added by aducharne 4 years ago.

Download all attachments as: .zip

Change History (7)

Changed 4 years ago by aducharne

Changed 4 years ago by aducharne

comment:1 Changed 4 years ago by aducharne

J'ai repris le code, fait des calculs, et l'analyse ci-dessus est cohérente avec le pb reporté par Arnaud Caubel, et mentionné dans les minutes de la réunion hebdo du 14/11/2017 : weekly meeting reports.

Des graphiques complémentaires sont dans le document http://forge.ipsl.jussieu.fr/orchidee/attachment/ticket/412/bug_ksat_Nov2017.pdf

Le bug est donc confirmé et entraîne une sous-estimation de ksat, d'un facteur 3.5 à 10 en surface pour les 3 classes de sol de Zobler, avec des conséquences potentiellement grandes sur tous les termes du bilan d'eau et la dynamique des débits.

La solution est simple, et ne concerne qu'hydrol_soil_coef, où le bloc ci-dessous devrait être remplacé par :

       DO jsl=1,nslm
          DO ji=1,kjpindex 
             
             ! it is impossible to consider a mc<mcr for the binning
             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji)))
             
             i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin)
             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), &
                  a_lin(i,jsl,njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji)))  ! in mm/d
          END DO 
       END DO

Seule la ligne portant sur k(ji,jsl) est modifiée, à faire aussi pour le bloc si ok_freeze_cwrr=T.
J'attache aussi le code hydrol modifié, à partir de [r4753].

comment:2 Changed 4 years ago by aducharne

J'avais oublié un point, c'est de corriger la sortie de la nouvelle variable ksat analysée par A. Caubel, et récemment introduite pour CMIP6 (r4545). Je joins en attachement de nouvelles sources pour hydrol.f90, qui termine je l'espère les corrections sur k et ksat, et qui supprime aussi le controel de la sortir de "kk_moy" par ok_freeze, pour avoir accès à cette variable importante (requise pour SP-MIP) dans tous les cas.

Le pb de fond, c'est que l'on ne calcule pas une bonne fois pour toute une valeur de ksat qui est modifiée par kfact et kfact_root, mais que ces effets doivent être intégrés à chaque fois qu'une utilise ksat ou k. C'est TRES buggogène. Ce sera à reprendre pour ORCHIDEE3.0 si on a le temps.

Changed 4 years ago by aducharne

comment:3 Changed 4 years ago by jgipsl

Done in trunk rev [4764]

comment:4 Changed 4 years ago by aducharne

  • Resolution set to fixed
  • Status changed from new to closed
Note: See TracTickets for help on using tickets.