Good afternoon,
I am using Hydrus 2D to simulate drywell infiltration in a 2D axysimmetric domain (x-z). The domain is made of 6 different horizontal layers whose hydraulic properties are known. The initial condition in the domain is imposed in terms of head, with each layer having a certain initial head.
By using the reservoir boundary condition, I am imposing a solute concentration of 10^(-5) g/cm^3 coming in the soil from the well. I am imposing a variable head at the well boundary, and the flow field does not give any issues (the whole thing works fine when I am not simulating solute transport).
The initial condition for the concentration in the domain is imported from the last time step of another simulation, and since the solute is supposed to represent oxygen, I have higher concentration at the top boundary that decreases with depth reaching its lowest value at the bottom boundary. Such initial concentration is in the order of 10^(-6) g/cm^3 (varying from ~8*10^(-6) to ~1.5*10^(-6)).
Oxygen can be replenished from the atmospheric boundary at the top, and in the domain I am imposing a zeroth-order decay rate for oxygen (representing microbial oxygen consumption) that is constant (and the same) for all layers.
The issue I am having is that despite injecting water with a solute concentration of 10^(-5) in a domain that has concentrations in the order of 10^(-6), numerical instabilities kick in so badly that I end up with concentrations of the order of 10^(-4) at the "interface" between injected water and "local" water (this starts happening from the very beginning).
I tried reducing the initial, maximum, and minimum time steps; I tried using different combinations of the space and time weighting schemes; I tried playing a bit with the stability criterion ( the product of Pe*Cr); I tried normalizing the concentration values (and the associated transport parameters) but I still and up with the same issue. The mesh is already very refined because of the "jumps" in hydraulic properties between layers. I am not sure what to try next.
For context, when the initial concentration in the domain is uniform an equal to zero, all other things being equal, the numerical instability is still present but not nearly as bad (the maximum concentration in that case is of the order of ~1.2*10^(-5)).
I would really appreciate any help with this, thank you!
Numerical instability for concentration
-
AlessandraBonazzi
- Posts: 5
- Joined: Fri Jun 07, 2024 8:44 pm
- Location: USA
Re: Numerical instability for concentration
I have to admit that I do not understand why you post this question here now after I have already looked at your project and provided you with the answer below:
Hi Alessandra,
I’m looking at your project. Here are some main comments:
1. The problem clearly is with zero-order degradation. If I remove this process, I can easily run your project (after a few minor modifications discussed below) within 4 minutes with mass balance errors below 1% for both water and solute.
2. Zero-order rate was always intended mainly for production, since you need to have an unlimited mass of some compound to sustain this reaction. With HYDRUS, it is usually used to simulate processes such as the mineralization of organic nitrogen (when simulating nitrogen fate and transport). As you likely understand, a zero-order reaction means that the rate is independent of the concentration of a given solute, i.e., it has the same rate whether the concentration is large or small. If you use it for degradation, the same rate applies whether the concentration is large, small, or even zero or negative, and therefore, you can clearly get negative values. Since this process was not intended in HYDRUS to be used for degradation, Hydrus does not stop this process when the solute is zero. Degradation is usually simulated as a first-order process, when the rate is proportional to concentration, and thus stops when the solute disappears.
3. If you insist on using a zero-order process, I could perhaps modify the code so that the process would stop once zero concentration is reached. However, I think you should first reconsider using this description of the degradation process.
4. You do not really use the reservoir BC. You use a time-variable head BC in the well and specify the water level in the well (with linear interpolation between specified values). If you want to use a reservoir BC, you will need to assign a seepage face BC in the well, specify the initial position of the water level in the well (which you may do), and then define water pumping or injection as a function of time in the well. I have posted a couple of examples on the HYDRUS website on how to use the reservoir BC.
5. You cannot have a free drainage BC on the side of the domain. Note that a free drainage BC implies a zero pressure head gradient and a unit gravitational gradient (i.e., q=K(h)), there is no unit gravitational gradient on the side and the flux is not q=K(h).
Other comments:
1. You certainly do not need to have such a fine spatial discretization. You should have fine discretization around the well, so that water flow solution converges. You can have coarser discretization further away from the well. The water flow solution should define the spatial discretization. Solute is usually secondary and should work once water works.
2. There is no need to limit the maximum time step. The maximum time step is determined by the model based on convergence of the water flow solution, and courant numbers for solute transport (to maintain stability of the solution). I would relax the max time step to 1 h.
3. I would limit the number of print times, unless you really need all these print times for some additional analysis. Note that your project is over 1 GB (likely the largest HYDRUS project I have ever seen).
4. Solute Transport- General Information: You should certainly consider tortuosity. It is a factor you cannot really neglect.
5. Solute Transport- General Information: You do not need to specify “Iteration Criteria” for solute transport, since solute transport is linear (and thus there is no need for an iterative solution).
6. I would increase ten times the dispersivity values. Note that you are dealing with a field process in undisturbed soils. Dispersivities smaller than 1 cm are common for repacked soil columns, but not for field soils. Additionally, your dispersivity values lead to Peclet number of 30, which may cause numerical oscillation (unless you use Upstream Weighting, which you do, but is likely not necessary).
7. Initial conditions: While water contents can be different in different soil layers and can change abruptly at material boundaries, pressure heads should be continuous. Pressure head is the driving force of water flow and, as a result, it needs to be continuous, even on material boundaries. It should not change abruptly (that’s unphysical). If you have measurements in different soil layers, then you should use these values at measurement depths, but then you should change the pressure head linearly between these measurement points. That would be way more physical than what you have now.
Hi Alessandra,
I’m looking at your project. Here are some main comments:
1. The problem clearly is with zero-order degradation. If I remove this process, I can easily run your project (after a few minor modifications discussed below) within 4 minutes with mass balance errors below 1% for both water and solute.
2. Zero-order rate was always intended mainly for production, since you need to have an unlimited mass of some compound to sustain this reaction. With HYDRUS, it is usually used to simulate processes such as the mineralization of organic nitrogen (when simulating nitrogen fate and transport). As you likely understand, a zero-order reaction means that the rate is independent of the concentration of a given solute, i.e., it has the same rate whether the concentration is large or small. If you use it for degradation, the same rate applies whether the concentration is large, small, or even zero or negative, and therefore, you can clearly get negative values. Since this process was not intended in HYDRUS to be used for degradation, Hydrus does not stop this process when the solute is zero. Degradation is usually simulated as a first-order process, when the rate is proportional to concentration, and thus stops when the solute disappears.
3. If you insist on using a zero-order process, I could perhaps modify the code so that the process would stop once zero concentration is reached. However, I think you should first reconsider using this description of the degradation process.
4. You do not really use the reservoir BC. You use a time-variable head BC in the well and specify the water level in the well (with linear interpolation between specified values). If you want to use a reservoir BC, you will need to assign a seepage face BC in the well, specify the initial position of the water level in the well (which you may do), and then define water pumping or injection as a function of time in the well. I have posted a couple of examples on the HYDRUS website on how to use the reservoir BC.
5. You cannot have a free drainage BC on the side of the domain. Note that a free drainage BC implies a zero pressure head gradient and a unit gravitational gradient (i.e., q=K(h)), there is no unit gravitational gradient on the side and the flux is not q=K(h).
Other comments:
1. You certainly do not need to have such a fine spatial discretization. You should have fine discretization around the well, so that water flow solution converges. You can have coarser discretization further away from the well. The water flow solution should define the spatial discretization. Solute is usually secondary and should work once water works.
2. There is no need to limit the maximum time step. The maximum time step is determined by the model based on convergence of the water flow solution, and courant numbers for solute transport (to maintain stability of the solution). I would relax the max time step to 1 h.
3. I would limit the number of print times, unless you really need all these print times for some additional analysis. Note that your project is over 1 GB (likely the largest HYDRUS project I have ever seen).
4. Solute Transport- General Information: You should certainly consider tortuosity. It is a factor you cannot really neglect.
5. Solute Transport- General Information: You do not need to specify “Iteration Criteria” for solute transport, since solute transport is linear (and thus there is no need for an iterative solution).
6. I would increase ten times the dispersivity values. Note that you are dealing with a field process in undisturbed soils. Dispersivities smaller than 1 cm are common for repacked soil columns, but not for field soils. Additionally, your dispersivity values lead to Peclet number of 30, which may cause numerical oscillation (unless you use Upstream Weighting, which you do, but is likely not necessary).
7. Initial conditions: While water contents can be different in different soil layers and can change abruptly at material boundaries, pressure heads should be continuous. Pressure head is the driving force of water flow and, as a result, it needs to be continuous, even on material boundaries. It should not change abruptly (that’s unphysical). If you have measurements in different soil layers, then you should use these values at measurement depths, but then you should change the pressure head linearly between these measurement points. That would be way more physical than what you have now.
-
AlessandraBonazzi
- Posts: 5
- Joined: Fri Jun 07, 2024 8:44 pm
- Location: USA
Re: Numerical instability for concentration
Thank you so much for the reply, I checked my email inbox and couldn't find your reply, I don't know what went wrong, but I'll work on my project following your suggestions.
Re: Numerical instability for concentration
I sent you my response on Fri, Sep 27, 2:24 PM (4 days ago), with cc to Scott and Veronica. J.
Re: Numerical instability for concentration
Bonsoir
lors de ma modelisation inverse, au niveau de Data for inverse solution, après avoir rempli le tableau j'ai eu ce message" the forth column of the record 1expo0004 must be smaller than or equal to number of observation nodes" alors que dans la colonne position j'ai mis les profondeurs de mes point d'observation. POUVEZ VOUS M'AIDER SVP A COMPRENDRE ?
lors de ma modelisation inverse, au niveau de Data for inverse solution, après avoir rempli le tableau j'ai eu ce message" the forth column of the record 1expo0004 must be smaller than or equal to number of observation nodes" alors que dans la colonne position j'ai mis les profondeurs de mes point d'observation. POUVEZ VOUS M'AIDER SVP A COMPRENDRE ?
Re: Numerical instability for concentration
Observation nodes are numbered from 1 to NObs (the number of observation nodes). It is this number (not the depth) that you have to place in the fourth column. J.