GROUP and Objectives

• Persons: Frederique, Agnes, Jan, Catherine, Fuxing, Tao, Philippe P., Jean Louis, Josefine
• Main responsibles: Frederique ?
• Objectives / roadmap : Possible roadmap (to be refined and detailed)
  • Integrate a « unique soil discretisation for W and E » (travail de Fuxing) => date? Current version of Fuxing is 8 Meters for both E and W, with same discretization. Questions: AD and JP wonder about the impact on the runoff and thus river discharges and possibly the time required for the spin up (such time is likely to increase with a 8 m soil) ?
  • Modify the coupling of the new multi-layer snow with soil to be fully implicit. => date ? The current version in MICT works (couple and force mode) but is not fully implicit and the energy conservation at each time step should be checked! The work to make it fully implicit needs to be quantified with Tao (Anticipated to be relatively straightforward)
  • Eventually implement a soil freezing for the "thin layers of the new vertical discretization" => ?? A solution to avoid large oscillations needs to be proposed and tested.
  • Improve the albedo of the vegetation and bare soil: Work that can be pushed by a Post-doc from the CMUG project (Patricia is responsible of the project; person to be recruited); => Aout 2015; propositions: Make some tests with the We can try to re-activate the dependance of bare soil albedo to sol humidy with the 11 layers scheme.
• Remarks : The agreement and the details of such roadmap should be done; Most likely we need to separate the different TASKs with a different group of people for each?

Discussions

Jan Polcher, 13/01/15 (mail)

après notre discussion cet après-midi j'ai un peu réfléchi à comment "merger" les 4 points de Philippe dans le trunk, tout en permettant à tout le monde de suivre ce qui se passe et donc de construire une confiance commune dans ces changements.

Je propose un plan en 4 étapes ou nous laissons quelques semaines entre chaque étape afin que tout le monde puisse faire des runs et des tests avec cette nouvelle version. Ceux ceux qui ont piloté les changement pourront écrire une documentation technique pendant ce temps. J'ai essayé aussi de séparer les changements structuraux des améliorations des paramétrisations actuelles. Donc voici mes 4 étapes :
1) implémentation de la discrétisation verticale commune entre T et W. L'important est de garder la profondeur en T et W comme des paramètres indépendants afin qu'on soit libre de choisir plus tard ce qui convient.
2) Enlever la neige de Thermosoil et la mettre en sandwich entre enerbil et thermosoil. Ceci est donc la partie thermodynamique de la neige que Thao vient d’implémenter dans MICT mais avec la l'interface neige/sol finalisée.
Après ces deux changements on aura une version qui a la nouvelle structure et les nouvelles approches numériques. Les paramètre auront peu changés et donc à priori le comportement du modèle sera similaire au trunk. Un environnement idéal pour voir si tout marche bien techniquement.
3) On rajoute les nouvelles capacité et diffusion thermiques de Fuxing avec le transport de la chaleur par l'eau.
4) L'albedo de la neige et les autres composantes du modèle de neige de MF sont rajouté afin d'arriver à la même configuration que dans MICT aujourd'hui.

Si on veut on peut rajouter une 5ème étape qui est l'activation du gel du sol et la vérification des processus dans la première couche.

Je pense qu'après ce processus on sera au stade ou voulait en venir Philippe pour CMIP6. Si à chaque étape on échange sur ce qui a été fait, que tout le monde fait les tests qu'il pense important on aura tous la même compréhension des évolutions faites.

Pour réaliser ce plan je suis disposé à travailler avec Fuxing et Thao sur les étapes 1 et 2. Je ne peux pas quantifier le travail pour chaque étape, mais une semaine me semble un bon first-guess. Je pense que pour la neige ce sera un peu plus long. Mais comme cela vient en second on peut préparer les choses en amont.

Philipe Peylin, 17/01/15 (mail)

Concernant la proposition de Jan:

• De mon coté je vois un interêt au point 1) et à le traiter tout de suite en discussion avec ce qui va devoir être fait pour inclure le permafrost (déja inclus dans MICT (version de Charlie Koven)) car dans ce cas la discretisation pour la Temperature devra aller beaucoup plus profond.
• Ensuite pour le point 2 et 4, personnellement j'aurais traiter cela ensemble (mais je ne suis pas le plus compétent pour juger de l'intéret de séparer)

Agnès Ducharne, 23/01/15 (mail)

My opinion is that the main challenges are :
1) to get the new snow scheme, as it improves the albedo, but in a way that fully respects the surface energy budget and the implicit coupling with LMDZ
2) to define the best vertical discretizations for T and W, in a way that keeps them as close as possible for convenience. The new soil termal properties are an improvement that can be added even if we do not change the vertical discretizations compared to the present Trunk, and heat convection might probably be overlooked in CMIP6.v1 if we cannot achieve identical discretizations, as it has a small effect based on Fuxing's results.
3) define "MODIS-like" reference albedos for the PFTs

I focus today on the second point: define the best vertical discretizations for T and W. Again, we can identify several issues, and from my point of view, the main ones are:
a) the deep soil required for T may be to deep to reasonable runoff values compared to observed river discharge
b) the very thin top layer required for W may be too thin for T and induce numerical instabilities of Ts and the surface energy fluxes
c) this behavior might be exacerbated in case of freezing

Attached is a series of plots by Fuxing illustrating point b ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/fuxing_np_point_numerical.pdf.

Explanations of the plots: If we look at the temperature at each time step, we can find oscillations over the upper 1.4cm soil depth, but the oscillations are not too large. Three simulations were compared (only 5-day simulation, at point: 25E,0N) (see attachment): CTL: the standard vertical discretization is used (2-meter for moisture, 5-meter for temperature), 7.5-min time step (black line) EXP: the revised vertical discretization (8-meter for both moisture and temperature), 7.5-min time step (red line) EXP2: same as EXP, but 15-min time step (green line) Main conclusions include: (1) in the revised soil layer (8-meter), the oscillations appear at upper 4 layers (0.5mm, 2mm, 6mm, 14mm). the oscillations disappear at 5th layer (29mm) (see Figs. c-g). But there is no very large or very small strange temperature in the sub-surface (Figs. a-b). (2) There are also oscillations in the FLUXLAT, FLUXSENS, SOILFLX, RAIN, QAIR (Figs. h-l), but the variations are in normal range. Besides, there are also oscillations by using the standard discretization.

The good thing is that there seems to be solutions to the above issues:
a) if pb a is proven (we are still waiting for an answer about this), we could do as in SURFEX-ISBA and only consider the same discretization over 2-m, below which T is solved under the assumption of a uniform W profile (consistent with teh gravitational drainage assumption). After talks with Jan, Frédérique and Philippe, I feel it's a solution anyone could agree on if cannot make deep W discretization work.
b) we know we have a pb here (Fuxing grahs), but Jan mentioned several times that numerical solutions exist. It's about filtering the variability, but I personnaly have no idea what is recover in practice, and it woud be nice if Jan could give us more information on this. Another solution might be to lump the 4 top layers of the W discretization to define the first layer for T It's seems rather easy, especially if we agree on eth proposition made by Aurelien during his PhD: stick to the present W discretization in the top 25cm, even we allow to change it below.
c) this still needs to be assessed, and we need to do it for CMIP6.v1 if we want to activate soil freezing.

Summary of a talk between Fuxing Wang and Aaron Boone at the AMA2015 in Toulouse

(1) Potential problem with the too deep soil depth: A. BOONE said that not only the annual mean discharge changes by using very deep soil layers, the seasonal cycles of discharge (phase and amplitude) are changed a lot by using very deep soil layers in ISBA model. And the second situation may be more serious that the first problem. It seems the hydrology model in MPI also has this problem mentioned by A. BOONE. Currently, the soil depth is 2-3m (changing for different regions) for hydrology,and it is 12m for thermal (constant) in ISBA model. The soil moistures below 2-3m are computed by Darcy's law.

Decharme, B., A. Boone, C. Delire, and J. Noilhan (2011), Local evaluation of the interaction between soil biosphere atmosphere soil multilayer diffusion scheme using four pedotransfer functions, J. Geophys. Res., 116, D20126, doi:10.1029/2011JD016002.

Decharme, B., E. Martin, and S. Faroux (2013), Reconciling soil thermal and hydrological lower boundary conditions in land surface models, J. Geophys. Res. Atmos., 118, 7819–7834, doi:10.1002/jgrd.50631.

(2) Spinup: For the spin up in the ISBA model, it is about 10 years (or more) for dry conditions, and it is much shorter for humid regions (2 year). And it is confirmed that the spinup length for temperature is shorter than that for moisture. For our revised model, A. BOONE commented that it maybe take several decades to get soil moisture equilibrium (but A. BOONE is also not sure exactly how many years are necessary).

(3) 1st layer depth It is 1cm for the 1st soil layer in ISBA. A. BOONE's opinion is that the 1mm soil depth is too thin, and A. BOONE's suggested it better to use at least 1 cm for the 1st layer.

Jan Polcher, 26/01/15 (mail) to Agnès 23/01/15

thank you for this analysis ... which I totally share. I will give you my feeling on a few of the points you raise :

The challenges :

• I have put the snow second because I find the interaction with Tao more difficult and thus I wanted to give us more time for it.

The depth : This should be configurable by the user as different applications will have different needs. We should code things so that each one can test the configuration which suits him/her best. The limiting conditions are :

• The top layers need to have the fine vertical resolution of the soil moisture.
• The total depth for soil moisture should be equal or smaller than for temperature (The total depth of soil moisture and temperature need to be parameters in the run.def).
• After that we can debate if some upper and lower bounds for the depth should be introduced.

As for the technicalities :

• soil moisture below its total depth will be take as constant and equal to that of the deepest layer.
• The smoothing of the soil temperature profile at depth of less than 5cm needs to be done in thermosoil once the new profile is calculated. That should be enough to avoid numerical instabilities.

For the snow, I am still waiting to see the equations to be used to coupled the bottom temperature of snow pack with the one of the soil. Once we have all the equations of the temperature diffusion (from the skin temperature to the bottom of the soil and through the snow pack) in implicit form, things should flow smoothly.

Fuxing Wang, 26/01/15 (mail)

(1) The depth. This should be configurable by the user as different applications will have different needs. We should code things so that each one can test the configuration which suits him/her best.

In the new discretization, the HYDROL_SOIL_DEPTH (=dpu_max), THERMOSOIL_NBLEV (=ngrnd) can be defined in run.def. I also added a flag (VERDIS_FLAG, can be defined in run.def ) for different discretizations: verdis_flag = 1, the standard discretization; verdis_flag = 2, the revised discretization. But 'nbdl' (Number of diagnostic levels) and 'nslm' (number of CWRR levels) in 'constantes_soil_var' (=11 in standard discretization) was not defined in run.def. Is this kind of configuration OK ? Do the variables 'nbdl' and 'nslm' also need to be defined in run.def ? Because they also change when using the revised discretization.

(2) The technicalities : The smoothing of the soil temperature profile at depth of less than 5cm needs to be done in thermosoil once the new profile is calculated.

In this case, we smooth the top 5 soil layers for soil temperature ? The bottom positions from top 4st to 6th layers are : 2.151cm, 4.497cm, 9.189cm.

Frédérique Cheruy, 27/01/15 (mail)

About the optimal vertical discretization for W and T
1- We still need to verify that the 8m deep soil, is a problem for the runoff and the drainage.
2-Afterwards, as already discussed last week, and reported in Agnes mail, we might want to have different depths for moisture and temperature but we want to have the dependance of the soil thermal properties with the soil texture (3 or 12?), and the same vertical discretization for moisture and temperature to avoid interpolations.
I am not sure we want to externalise nslm or ngrnd, and let the user play with it. If it is an advanced user he can go into the code and change the value of the parameter, if not may be it's better he does not play with it at all.
But we need to test upper and lower bound that we want to adopt for CMIP6 at least.
Once we get the answer to point 1, Is the 8m T, 2m W with same vertical discretization for the upper 2 m and moisture constant and equal to the one of the deepest layer a good point to start with? Can we agree on that?
3- Is someone able to explain to me (sorry ), why a smoothing is a good way to handle with numerical instabilities, to me it looks like an artificial way of masking a potential problem.
4- Do we want to consider the problem of the soil freezing now?

About the albedo
Whoever will work on that subject, it would be good to know if the bare soil albedo maps prepared by Sebastian and team can be directly pluged into Orchidee-trunk or if we need to see with ISBA people, or to redo the work ourselves.

Jan Polcher, 27/01/15 (mail) to F. Chéruy, 27/01/15

Je suis d'accord avec toi qu'il faudra définir la "configurabilité" des niveaux verticaux en fonction des différentes contraintes. Le nombre de nivaux sera sans doute un résultat plutôt qu'un paramètre.

Sinon pour le filtre. Tu as une onde thermique qui a une vitesse de propagation dans le sol. Il faut avec la discrétisation satisfaire le critère de Courrant-Levy : ne pas sauter plus d'une maille par pas de temps. Cela n'est pas possible avec les niveau de l'hydrologie. Au lieu de filtrer les ondes rapides (On ne veut pas de cela !) il faudra filtrer le résultat et aboutir à une discrétisation effective (au plan numérique) plus lâche que celle de l'hydrologie.

C'est la même chose qui est fait dans LMDZ au pôle !

Jan Polcher, 06/02/15 (mail)

I have given it some thought on how to specify the parameter in order to choose the common levels between hydrol and thermosoil.

In order to see what works best from a user point of view, I have written a simple python script to select the properties of the vertical discretisation and plot the result. The 4 parameters I propose for the moment are :
Depth of top hydrology layer [9.77517107e-04 m] >
Depth of hydrology [2.0m] >
Depth at which linear hydrology layers start [0.75 m] >
Depth of thermodynamics [8.0m] >
in [] are the default values you can if you do not enter anything to the script and which produce a result very similar to what we have today. Attached is also the plot produced by the script ​http://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/VerticalLayers.pdf. The script which allows to plot the desired configuration is given here ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/NewHydroVert.py.

There is still the open question if we need a vertical filter for the top layer of thermosoil, and if so we need to specify the filtering depth.

So please play with this and see if these are the parameters which are best suited for ORCHIDEE or if we should change some.

The way I have coded the calculation of the vertical levels and number of layer in Python, are very FORTRAN like. So once we agree it should be easy to put that into a module of ORCHIDEE and test again if we need or not a numerical filter.

So once we agree on the choice it can be implemented in the trunk.

Agnès Ducharne, Fuxing Wang, Fréderique Chéruy, 06/02/15 (meeting)

We had a skype meeting today about the soil discretization for water and T. We looked at three discretizations (D1: as in Trunk; D2: 8m for water and T, with the same 17 nodes, and the same location of the top 11 nodes as in D1 for water; D3: 2m for water as in D1 + assumption of uniform profile below, 8m for T as in D2, so that the two diffusion schemes use the same nodes in the top 2m)

• Spin-up: Fuxing tested this in off-line mode, by repeating the same year, globally (forcing dataset: ergon/IDRIS: /linkhome/rech/psl/rpsl035/IGCM/STORAGE/BC/OL2/CRU-NCEP/v5.3/twodeg/cruncep_twodeg_1990). Almost all land points reach equilibrium in 20 years for D1, and you need more than 40 years with D2. Fuxing still needs to check how long it takes to reach thermal equilibrium in D3, then to do the same tests online.
• The goal of the spin-up phase is to compare the different simulations in a comparable state: the next steps will be to compare the water budget + seasonal cycles of T, top SM, surface fluxes and river discharge, between D1 and D2, and the thermal state between D1, D2 and D3. To be tested off-line then on-line. This should allow us to decide which kind of vertical discretization is preferable for CMIP6.v1. No conflict with the above suggestion by Jan, which makes the vertical discretization more flexible.
• Then, we'll have to check how the chosen discretization performs with the developments of Fuxing regarding the soil termal properties (dependence on soil texture) and heat convection (from which we don't expect a large effect).
• We shall not forget the issue of the large T variations in the top soil layer in D2 and D3, with a thin top soil layer for heat diffusion.
• If the group agrees on using the Reynolds map of texture for CMIP6, with 12 texture classes, this will impose to have a soil depth of at least 9 m for T (then why not 10 m?) + basic checks on the water budgets and T behavior given the new textures.
• Shall we include permafrost for CMIP6.v1 ? with which soil depth and discretization for T ?

Frederique Cheruy after Skype meeting and discussion with Jan ( 7/02/2015)

• Spin-up and off-line simulations.

In FG2low off-line simulations, the SM patterns appears to be very strange. We suspect it is due to the forcing itself. It would be very helpful to have off-line simulations with an other forcing.

Can someone recommend a forcing data set and make it available for the spin-up tests with the new vertical discretization?

Fuxing Wang, Fréderique Chéruy, Agnès Ducharne, Jan Polcher, Philippe Peylin, Jean-Louis DUFRESNE, et al. 11/02/2015 (mail)

Updated history: 17/02/2015

Topic:

Oscillations in the ORCHIDEE-LMDZ coupled model with the revised soil thermodynamics.

Background:

There are mainly three modifications in the soil thermal process.
(1) In the current ORC, different soil vertical layers are used for hydrology (11 layers, 2-meter) and thermal (7 layers, around 5-meter). In the revised one, the same vertical discretization (8-meter) is implemented for both models;
(2) In the current ORC, the soil thermal conductivity and soil heat capacity are not changing with soil textures. In the revised one, the soil texture effects are included;
(3) the soil heat convection (by liquid water transport) is neglected in current ORC. The soil heat convection and conduction are coupled in the revised ORC.

Models:

The LMDZ5B and ORCHIDEE trunk are used. The time step is 7.5 mns, but the radiation is called at each 1 hour time step.

Experiments:

See Table 1: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/Table1_oscillations_LMDZOR.pdf

Results:

Two points under different climate conditions are selected: one point at (25E,0N), humid; another point at Sahara (25E,15N), arid.
Please see (together with Table 1):
Figs. 1-6 : ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_general.pdf,
Figs. 7-10 : ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_radia_step.pdf
Figs. 11-12: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_radia_step_atm_25e0n.pdf
Figs. 13-16: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_radia_step2.pdf
Figs. 17-18: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_cre2.pdf
Figs. 19-22: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_cdrag.pdf
Figs. 23-24: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_nbapp192_lwcpl.pdf
Figs. 25-26: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_lwdown.pdf
Fig. 27: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/test_oscillations_orc_revised_thermo_cre.pdf

Conclusions:

(1) When the radiation is called every hour, the oscillations originates from the CDRAG and the coupling with the atmospheric PBL.
(2) When the radiation is called every time step, the oscillations that may appear originate from the cloud and the PBL schema.
(3) The lwup has almost no impacts on the oscillations; the oscillations are alleviated with lwdown (new). but they are not removed.
See Table 1: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/Table1_oscillations_LMDZOR.pdf

Frederique Cheruy Summary (11/02/2015)

I am trying to summarize where we are now:

1- The numerical oscillations seem to be more on humid area. 2- They appear to be pulsed with the update of the radiation. This is true for the old vertical discretization as well as for the updated ones. However the amplitude of the oscillations are larger with the updated vertical resolution. 3-When the radiation is updated at the same frequency as the physics (7.5 mn),

It is worse (whatever the vertical resolution for the soil-thermics is), which makes me to state that

the probleme is not due to the delay in the update of the radiation.

The SW net radiation has numerical oscillation as well, which are related noisy cloud cover and noisy convection related variables.

My feeling is that the problem exists independently of the soil vertical discretization, but it is somehow amplified when the upper soil layers are thinner.

1- I would suggest to look in detail at the atmospheric tendencies (temperature), for convection, thermics ., verical diffusion ... 2- However, I do not understand why calling the radiation with the same frequency as the rest of the physics makes things worst.

Has someone an other suggestion? Frederique

Mail JLD

Hi Frederic,

I spend some (short) time with Fuxing this morning, and we have the same filling that the origin of the oscillations is not in the soil but in the boundary layer and the clouds (and they are amplified when the surface reacts quickly). Fuxing is doing some complementary tests.

J-Louis

Mail Pascal Maugis: functionning of CWRR

Dear all, thank you for associating me to the discussion. I paste here what I wrote to Agnès yesterday on the discretization scheme of the CWRR hydrology. Note that this comes exclusively from Patricia de Rosnay's thesis and that I did not check whether actual implementation follows scrupulously the given matrix equation (although it would not be very long to check).

Il y a donc une relation : la loi de Darcy, écrite en Différences Finies du premier ordre, avec pour variables les flux Q et la teneur ponctuelle Thêta aux noeuds ; et une équation : la conservation de la masse, avec pour variable la quantité d'eau W comprise entre les deux demi-intervalles autour du noeud (la maille), non affectée à un point particulier mais référencée avec l'indice correspondant au noeud. W peut aussi être vu comme l'intégrale d'une teneur en eau locale continue et linéaire par morceaux entre les noeuds avec pour valeur les Thêtas aux noeuds. La teneur en eau moyenne sur la maille ne coïncide donc pas avec la valeur au noeud. Noeud qui n'est d'ailleurs pas le centre de la maille. Avec une progression géométrique des intervalles d'un facteur 2, cela place le noeud au premier tiers de la maille. Les rapport 3/4 et 1/4 qui ont peut-être été mentionnés ne correspondent pas à des positions de noeuds, mais à des poids entre les valeurs de Thêta pour obtenir la moyenne de la teneur en eau sur chaque demi-intervalle. Donc si tu veux représenter la teneur en eau, ce serait numériquement cohérent avec l'interpolation de la formulation numérique que de relier les valeurs aux noeuds par des segments de droites.

Comme autre conséquence de ce choix de solution mathématique (l'espace des fonctions 1D continues et linéaires par morceaux), le flux est constant entre les noeuds (DF classique), et discontinu aux noeuds où le bilan est écrit et réparti sur la maille, et décrit par la variable W et ses termes sources.

On a donc des noeuds, quelque part dans les mailles, et des mailles délimitées par les moitiés des intervalles adjacents aux noeuds. On peut voir les choses alternativement en disant que nous avons des mailles séparées par des interfaces et que les gradients sont calculés par les différences de valeur de la teneur en eau moyenne affectées en un point particulier des mailles (avec une exception pour les mailles aux limites) : ça ressemble très fortement à du VF avec la méthode diamant (qui elle prend la valeur aux centres des intervalles). D'où une certaine confusion dans la description du schéma, qui est bien du DF simple mais qui peut être aussi interprété comme du VF diamant particulier.

In english in a few words, the formulation seems to be Finite Differences (FD) of first order. Nodal soil wetness Theta indeed corresponds to the numerical interpolation, which is linear between node (so I understood). Meshes apply for water balance, with variable W, and are delimited by half-way points between nodes. So the interpolation of Theta is not linear inside the mesh, since the node is inside it (at 1/3 of its height starting from the top). Water flux Q is constant between nodes, and their difference beween top and bottom of the mesh yields the water balance (storage + source/sink).

From what I understand of your intentions, I can say general things, depending on the way the coupling is considered. To avoid numerical inconsistencies stemming from using variables infeodated to the numerical formulation of one model (hydrology / thermodynamics), the coupling should rely on state variables and non-linear parameters only. Cross-use of fluxes should be avoided. There are TWO discretization, and therefore, TWO different interpolations to make between Water and Temp. variables:

• variables used to calculate the parameters related to water (resp. energy) balance that apply to the whole mesh upon which they are assumed constant. For water, it's only those used to calculate Tr_i, the integral of the sink term on the mesh.
• variables used to calculate the fluxes between nodes (assumed constant in FD). For example, if D depends on Temperature T finely discretized otherwise, you choose whatever relation to produce a single value of the equivalent D that will apply between Theta nodes (I will call it a water layer, not the mesh). You can choose the average of D over the water layer, or the value at the top node, or bottom node, or central node, etc. The difference will be in the dynamic of the coupling, which is linked to the intensity of the feedback and the general direction of the coupling. For example, if you want more dynamic of the downward thermal flux on water flux, use the top value. For more versatile use (which I imagine is your case since freezing and thawing come from top and bottom at different times), I would suggest getting an internal value : either the average of D(T) over the water layer, or the value of T interpolated at the center of the mesh according to the thermal formulation, or - less regular - the value at the water node.

Same thing could be done for use of parameters of the thermal equation influenced by water state variables.

Now this discretization issue is not necessarily connected to the oscillation you notice. It can also be due to the non-linearity of the relationships yielding the parameters of the equations ; be them internal to an equation, or coupling across equations. If a parameter (let's say the thermal diffusivity lambda) is highly sensitive to a state variables (be it T or Theta), small variations of T or Theta can produce large variations of T. Then, the only way is to buffer those relationships, either spatially by average over several cells (as it has been proposed), either temporally using sub-relaxation (only a fraction of the new value is used to update the parameter), either by modifying the relationship itself by smoothing it, or all of them ! Changing the unknown is also practiced in hydrology (switching to capilary pressure at high saturation) as well as introducing buffer terms on both sides of the equation. Actually, on such a numerical issue, imagination and trial-and-errors prevail. I have also to mention what I was conducted to do - like many other numerical modellers by the way - : lowering the time step. So the jumps in parameter values are lower. And do multiple iterations of coupling. However, this is costly (except with advanced high-time-order methods).

As you know, the key to understanding how the oscillations appear is to track which single parameter begins first to change dramatically, and how it propagates to the unknown, inducing oscillations maybe a few time-steps later.

I hope these few considerations will be useful.

Sincerely,

Pascal.

Note on soil freezing by Isabelle G. (email exchange with Philippe P)

Non, je n'ai pas désactivé le gel pour les petites couches. En revanche, j'ai créé une variable qui stocke en mémoire l'énergie libérée au cours du gel, sur l'ensemble du profil. Au moment du dégel, une fois que le sol a complètement dégelé, on compare ce stock d'énergie à l'énergie consommée pour dégeler le sol. Si cette énergie est supérieure à la première, on considère qu'on a trop "fondu" et qu'on aurait dû réchauffer à la place, donc on réajuste la température, sur l'ensemble du profil.

Ta prochaine question sera sans doute pourquoi j'ai travaillé sur l'ensemble du profil et pas par couches. La réponse est que j'avais la vieille version de la neige, implicite dans le sol. Donc, les "couches de sol" changeaient dans le temps en fonction de la hauteur de neige, donc il était difficile de garder une trace de l'énergie du gel/dégel du sol par couche.

Isabelle

Meeting on discretization and snow coupling issues (26 Fev 2015)

Persons: Catherine, Agnes, Frederique, Joséfine, Jean Louis, Gerhard, Tao, Fuxing, Jan, Philippe

Location: Jussieu

Vertical Soil Discretization

Presentation (Fuxing)

Fuxing described his work on the soil W and T discretization: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/20150226_Fuxing.pdf

He tested:

• standard discretization : V1
• new common scheme with 8m for both W and T: V2
• new common scheme with 2m for W and 8m for T: V3

Some results for T:

• 5 meters is not enough for diunal/seasonal cycle of T (criterial is that the amplitude of the cycle at lowest layer should lower by a factor "E-3" of that of the first layer ???)
• we need a depth of at least 9.5 m for sand soil when using new USDA texture classes.

For soil moisture:

• Overall the differences between V2 and V3 are not too large for the fluxes ?
• V2 need a longer spin-up (at least 40 yr compared to 20 yr for V3). Soil moisture at t0 = 0.2 m3/m3
• Impact of the drainage is large for JJA, but additional graphics are needed to better quantify the changes
• The impact on the river dischages is only significant for some rivers: additional test ?

Oscillation of at the surface:

• Fuxing did not had time to present the results but comprehensive results in his presentation and above discussion
• Most likely the origin of the oscillation is independant of the surface scheme and comes from the atmospheric PBL;
• The smaller the soil surface T depth the larger the oscillation seems to be on average
• More work will be done by Fuxing investigating causes from the PBL

Discussion

In general the impact of discretization on Temperature is small

The impact of the depth on Water balance seems less important than in ISBA but more checks are needed

We discussed how we should interpret the "W scheme" and the differences between the "Theta" values at the nodes (that are not the center of the layers) and the "water-content" values that are the integral of the water content for a given layer (and that are thus associated to any point). The specific case of the surface layer with the first node at the surface needs to be keep in mind.

Link between W and T scheme: for the heat capacity we should use the mean water content of the layer; for the conductivity, we should preferably use the water content at the interface. However, any choice is not that critical given the uncertainty on the physical laws.

Soil Freezing:

• change significantly soil heat capacity;
• we cannot anticipate the energy needed to freeze water as we dont know at the begining of the time step if the temperature will cross the zero before calculating it (thus we can not anticipate the release of energy).
• the scheme of Isabelle Goutevin has been implemented in the TRUNC.
• As for energy conservation the scheme buffers the energy released during freezing for all layers, then calculates the consumed energy during the next thawing period. If the the two energy do not match she distribute the excess/deficit of energy by re-adjusting the Temperatures of the whole soil profile. There is thus only energy conservation throughout the season (see email above by Isabelle). She could not do it per layers as with the old snow scheme the depth of the layers varies in time.
• With the new snow scheme we can envisage a new scheme that could conserve the energy each day or time step and possibly for each layer. To be discussed; but we will first use the scheme of IG as it is (for CMIP6).

Proposed Scheme

The objective is to propose a scheme flexible enough to allow users to test different configurations depending on their objectives and to be able to account for deep soil thermic representation in the case of permafrost.

Scheme (to be improved by all...): ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/soil_discretization.pptx

The main principles are:

• The T scheme is always as depth or deeper than the W scheme.
• We plan already a scheme that can account for large depth needed for permafrost in particular (more than 50 m)
• W and T share the same common discretization for the common part: the layers are the same; For the nodes: W has the firt node at the surface while T has it at the middle of the first layer; Else they have the same node position; Exception may also concern the bottom of W scheme: should we use Ti at the same level (bottom of the layer) ?
• The upper part of the scheme follow a geometric discretization of reason 2 (like the current W scheme)
• Below a certain depth, depth_lin, we have for the W and T a linear discretization, with the depth of each layers equal to that of the lowest layer of the "geometrical part" of the scheme. Note that depth_lin may be equal to depth max with thus no linear discretization section.
• Below another depth, depth_geom, we have again the geometrical discretization with a reason xx to be defined. This allows to model a deep soil for permafrost modelling. Note that Depth_geom may also be equal to depth_max with thus no bottom geometrical discretization.
• For the NUTRIENTS cycles, we will adopt a principle: always gather the existing layers to define a new vertical discretization adapted for C/N/P cycles.

The list of parameters to define the scheme would be:

• depth_max : maximum soil depth for W and T
• depth_lin : depth below which linear discretization occurs
• depth_geom : depth below which geometrical discretization occurs again
• reason_geom_up: equal currently to 2 (but we could possibly change it)
• reason_geom_below: we should decide what is the best following the test for permafrost
• depth_Wmax : max depth of the W discretization, between [depth_lin, depth_geom]

Diagnostic variable of the discretization would be:

• the number of layer for W and T scheme : Name to be defined ?
• the depth of the nodes for W and T
• the depth of the layer limits
• ??

Modulation of the Thermal properties as a function of the Water content:

• Heat capacity of each layer (C) should be a function of the layer water content: thus a function of the layer total water content ("W" variable of CWRR divided by the layer depth)
• Conductivity between layers (Lambda) should be a function of the "Theta" of layer i and i+1 (weighted arythmetic mean using the depth of layers)
• For the part of the thermal scheme that is below the hydrology scheme: we should take a constant in time value of "Theta" and "water content", to avoid changing heat capacity over time and thus not conserving the energy (without any corresponding imput/output of energy by water flow); the values to choose is open for discussion/tests ?

Note that proposed scheme includes what Fuxing has implemented so far and allow to account for deep permafrost layering as done in MICT

Snow coupling

Presentation by Catherine and Tao (MICT version)

We had time only for a short presentation of the new snow module implementation: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/Orchidee_snow_CMIP6.pdf

Some features :

• the number of layers of snow module (ES) is always 3
• when there is a snow fall the whole grid is covered by snow ; only the albedo is calculated with a fractional snow cover (standard linear fractional cover)
• After "Enerbil" the "ES" module update the snow temp profile and accounts for melt, freezing, sublimation,... It then computes the coeficient for the next time step that is given to "Enerbil": the equivalent capacity "C" term that account in an implicit way the ground heat flux (here the flux through the snow pack); Then he provides "Thermosoil" with the lowest snow layer temperature. Thermosoil use that temperature to update its temperature profile and provide back a surface flux for "ES" for the next time step.
• It seems thus that the scheme is not fully implicit with respect to the "flux from ground to the snow" that is provided in an explicit way from the previous time step value.
• With the coupled simulation, Tao had to define a minimum snow depth (1 cm ?) to avoid "crashes"; possible problem of energy conservation ?

Discussion

We discussed the potential impact of the non fully implicit "soil-snow-atm" scheme. Catherine mentionned that in ISBA, Aaron B. had similar issues and that the change to fully implicit did not impacted much the results.

However, we agreed that with very thin snow layers this could be a problem (see "trick" of Tao in coupled mode).

We proposed that a first improvement is to account for fractional snow cover for the energy budget (not only for the albedo). For that the equavalent "C" term passed to "Enerbil" should be a mix between the standard one of "Thermosoil" and the one from "ES". But this needs some revision of the equation and calls of the routine (to be writen explicitely first).

JP propose that we also split the fluxes from Thermosoil back to enerbil and ES with different values.. this would ensure also a full implicit scheme.

PROPOSED ACTIONS

• Vertical discretization:
  • Fuxing will correct the new scheme to account for the remarks on the heat capacity (C=f(wi)) and thermal conductivity ("landa" f(theta_i, theta_i+1))
  • Fuxing will provide additional diagnostic for the soil water balance (runoff, drainage, ...) to check the differences between 8m and 2m W depth; BUT for that Fuxing NEEDS a better forcing than the current one (ongoing action from the whole ORC group)
  • Jan will propose a little python code that visualize all possible discretizations: for us to check
  • we all review/agree/comment (for a restricted 2-3 weeks time period) on the new proposed discretization
  • FOR CMIP6: we propose to start with version V3 of Fuxing unless new tests show that V2 is better; we also proposed to use new USDA soil texture classes.
  • Remark on the W depth: in CMIP5 the depth of choisnel was increase to avoid Amazon drying; this should not be the case with the 2m CWRR as there are much more water holding capacity than with Choisnel. BUT we keep in mind the possibility to increase Depth_Wmax
• Snow Scheme:
  • Josefine with Tao will implement the updated version of ES (currently in MICT) into TRUNC
  • We then further add a fractional snow cover so that the problem of "very thin" snow layer (linked to the non fully implicit scheme) becomes most likely negligeable. FOR THAT, we need to start an email discussion about the requested changes (how does it change the call sequences to "thermosoil" and "ES" ?): Anyone is wellcome to launch the email discussion. Then, JLD propose that Fuxing helps with this implementation but other organizations are possible.
  • Few persons needs to gather to establish properly the equation-changes that are needed for a full implicit implementation between "enerbyl" "ES" and "Thermosoil": the list of the "few" needs to be decided ? (volunter ?)
• Soil freezing
  • We start for CMIP6 with the current scheme of Isabelle G.
  • GK/PP will ask IG if she wants to further work on improving the scheme (given that she is working a bit with ORC) with the objective of energy conservation per layer and not only over the whole year; Other person to work on it are welcome!

Comments on the vertical discretisation by Jan

I worked on the vertical discretization scheme as we outlined it in the meeting. I coded up the scheme we chose and found a few problems which we oversaw during the discussion.

The list of parameters as I see them today :

• depth_max = OK
• depth_Wmax = OK
• depth_geom = OK
• reason_geom_below = OK

• depth_lin : is probably a bad choice of words. We decided it to have a region of layers with a constant thickness. So I propose the following name instead : depth_cstthickness.

• depth_topthickness : this is a parameter we forgot. In order to start the geometric series, without a fixed number of layers, we need to choose the thickness of the top most layer. Furthermore this is an important parameter as it is key in the infiltration of water in CWRR.

• depth_refinebottom (Boolean) : Aurélien and Tristan have demonstrated that if we want to change the lower boundary condition in CWRR we will need to refine the resolution at the bottom. Using symmetry, we can easily set a refinement at the bottom while respecting all the above parameters.

There are a few dependences between these parameters. I have tried to put them into the code as well.

This is all coded up in python for testing purposes( ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/FinalHydroVert.py). It prints to the screen the levels computed and makes a graphic of the levels for T and W and corresponding layer thickness. So please play around with it and if all are happy we can implement it in the code.

Revised version of the vertical layering

After a meeting between Agnès, Frédérique, Fuxing and Jan (6/05/2015), the attached formulation of the vertical layering was agreed to. This is an implementation of what we understood from the CMIP-6 meeting.

Please try this code and provide comment so that it can be integrated into the trunk : ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/Revised_HydroVert.py

A decision was made to apply to the depth of the last hydrological node the same definition as the top one. That is, it is considered to be place at the bottom (or top) of the layer for numerical reasons. Thus the temperature node is put in the middle of the to and bottom hydrological layers.

Testing different soil vertical discretizations

26-05-2015.

Model:

ORCHIDEE Trunk (no STOMATE, with river routing)

Forcing:

CRU-NCEP/v5.3.2/twodeg

Definition of vertical discretizations (VDs):

VD1: 2m for water, 5m for temperature (trunk)
VD2: 8m for water and temperature
VD3: 2m for water, 8m for temperature

Simulations:

climatology simulations by repeating 1980-1989 for several loops
VD1: 5*10 years
VD2: 9*10 years
VD3: 5*10 years

Original data on dods:

VD1: ​http://dods.idris.fr/rlmd981/OL2/PROD/coseth20/COSETH20/
VD2: ​http://dods.idris.fr/rlmd981/OL2/PROD/coseth14/COSETH14/
VD3: ​http://dods.idris.fr/rlmd981/OL2/PROD/coseth30/COSETH30/
(in VD3, the forcing between different decades has some differences because the interpolation of different calendar between CRU-NCEP and ORC)

For those who have an account on ERGON/IDRIS, the data can also be found from here:
VD1: /linkhome/rech/lmd/rlmd981/IGCM_OUT/OL2/PROD/coseth20/COSETH20/SRF/Output/MO
VD2: /linkhome/rech/lmd/rlmd981/IGCM_OUT/OL2/PROD/coseth14/COSETH14/SRF/Output/MO
VD3: /linkhome/rech/lmd/rlmd981/IGCM_OUT/OL2/PROD/coseth30/COSETH30/SRF/Output/MO

Results:

(1) Checking whether the spinup criteria (for soil moisture SM and soil temperature T) has been satisfied with different VDs.

​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/spinup_histograms.pdf​

Critetia1 for SM: (TMCend - TMCstard)/TMCmean, with ±0.05 as threshold
Critetia2 for SM: (MCend-MCstart)/0.005, with ±1 as threshold.
Critetia for T: (Tend-Tstart), with ±0.1 as threshold.

(2) The comparison of E, D, R, Ts, and H for the 3 VDs after spin-up (Annual, January, July, Amplitudes)

(E: evapotranspiration, D: drainage, R: surface runoff, Ts: surface temperature, H: sensible heat flux)

Annual: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_SurfMet_Annual.pdf​
January: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_SurfMet_January.pdf​
July: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_SurfMet_July.pdf​
Amplitudes: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_SurfMet_amplitude.pdf​

(3) The evaluation of evapotranspiration with other datasets (climatology)

​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_Eobs.pdf

The ET products include Jung, LandFluxEVAL, Fisher, and GLEAM.
The maps are plotted for different seasons (MAM, JJA, SON, DJF, Annual).
Figs. a-e are the ORCHIDEE reference value.
Figs. f-j: ORC - Jung
Figs. k-o: ORC - LandFluxEVAL
Figs. p-t: ORC - Fisher
Figs. u-y: ORC - GLEAM
Zonal means at the 4th page !

(4) The comparison (climatology) of hydrographs between ORCHIDEE (VD1, VD2, VD3) and GRDC data

ORCHIDEE data: the last decades of spin up for each VDs.
GRDC data: Variable 'mergedhydro' in 'GRDC_Monthly_June14.nc'.
Selected basins: the catchment area > 1.5e6 km2.

Two different scripts are used to evaluate the hydrographs:
(a). A script written by Fuxing.
The calculation method:
First, identify the coordinates of the GRDC stations with a catchment area larger than 1.5e6 km2 .
Second, find the hydrographs of ORCHIDEE at the identified coordinates and plot them
​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_hydrographs_grdc.pdf​

(b) A more intelligent Python program developed by Jan.
​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_hydrographs_grdc_Jan.pdf​
and the river description: ​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_RiverDescription.txt​​

The station order in the figs of above (a) and (b) is different, but the line colour and the figure title format are the same for each station.
Notes: The hydrographs calculated from the two different program (a and b) are close:
--> the hydrographs of ORCHIDEE is underestimated over most regions, except for 'Parana do Careiro' which is overestimated.
--> it is significantly underestimated over OB, Mackenzie, Lena, Yenisey.
--> it looks fine over Congo and Rio Pahana, but still large discrepancy (underestimated) over Mississippi, Changjiang, and Amazona.
--> the hydrographs of VD1 and VD3 are close, but not for VD2.

Another simulation did by Sarah:
It is with the WFDEI corrected with GPCC data forcing file (at 3 hourly time step)
​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/Sarah_debits-ob-yenesi-lena.pdf​

(5) The snow mass by ORCHIDEE

Comparing over the last decades of each VDs.
​https://forge.ipsl.jussieu.fr/orchidee/attachment/wiki/Meetings/CMIP6/Physic/vd1_vd2_vd3_swe.pdf​​
Notes: Sorry I made a mistake on the SWE unit when plotting the original maps.
It has been corrected. Now the unit of SWE is the same with that of ORCHIDEE standard output (Kg/m²).