| | 263 | == Mail Pascal Maugis: functionning of CWRR == |
| | 264 | |
| | 265 | Dear all, |
| | 266 | 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). |
| | 267 | |
| | 268 | 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 ; |
| | 269 | 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. |
| | 270 | 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. |
| | 271 | 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. |
| | 272 | |
| | 273 | 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. |
| | 274 | |
| | 275 | 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. |
| | 276 | |
| | 277 | 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). |
| | 278 | |
| | 279 | 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. |
| | 280 | There are TWO discretization, and therefore, TWO different interpolations to make between Water and Temp. variables: |
| | 281 | - 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. |
| | 282 | - 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. |
| | 283 | Same thing could be done for use of parameters of the thermal equation influenced by water state variables. |
| | 284 | |
| | 285 | 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. |
| | 286 | 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). |
| | 287 | |
| | 288 | 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. |
| | 289 | |
| | 290 | I hope these few considerations will be useful. |
| | 291 | |
| | 292 | Sincerely, |
| | 293 | |
| | 294 | Pascal. |
| | 295 | |
| | 296 | |
