Recharge via time-dependent Neumann & deliverability boundaries in EOS9
Hello,
I am trying to model recharge to a shallow aquifer using a 20 m-deep 1-D (1 m x 1m cross-section) iTOUGH2-EOS9 model in which I incorporate the following processes: (1) infiltration at the top of the unsaturated zone, and (2) flow out of the bottom to represent saturated flow away from the model location. I initialized the system with the water table at 6 m deep by setting a hydrostatic water pressure profile below that elevation and gravity-capillary equilibrium above it.
I then ran a "steady state" simulation in which I incorporated a constant Neumann boundary in the top node to represent a yearly average infiltration rate, and a deliverability boundary in the bottom node to represent saturated flow. For the bottom boundary, I set the flowing bottomhole pressure to that of the bottom node from the hydrostatic initial condition (so that flow would only occur when additional water reaches the saturated zone and thus increases the pressure at the bottom), and adjusted the Productivity Index (PI) until the fluxes through the two boundaries were somewhat consistent. This resulted in what seems to be a rather low value for PI (on the order of 1e-17).
I am now having trouble when I try to incorporate transient infiltration using a simple "step function" boundary at the top, with a non-zero flux for some time and then zero flux thereafter. The simulation runs, but the results don't make sense because neither the total liquid mass in the system, nor the saturation profile in the unsaturated zone, nor the water table elevation, change through time. Could someone please help me to solve this problem by looking at the 3 attached files and answering the following questions?
- Why doesn't the code honor the flux-versus-time information I give it? According to the GENER block in my input file called "run.txt", the flux at node A0000 should be 0 at t=3601 seconds, but the output file called "run_gener.csv" shows that the rate becomes half of its initial value at 7200 s, and doesn't reach 0 until 10800 s. Why isn't the code taking more than one time step between each output time that I list in the TIMES block?
- Was the method I used to set the Productivity Index reasonable? The value I reached is ~13 orders of magnitude smaller than the values given in some of the examples in the TOUGH manual, but since I'm not actually modeling a well, I'm not sure what a reasonable value might be. Likewise, does it make sense to apply this boundary to only the bottom node? If not, how does the code interpret the input for the "number of layers"? That is, would I need to set the boundary at the top of the producing zone, or the bottom?
Thank you very much in advance for your help.
Mike
5 replies
-
Mike,
For #1, check out MOP(12), which you appear to have set to 0, which is a triple interpolation scheme. It sounds like you would like something closer to MOP(12)=2, which is the rigorous interpretation.
Let me know if that helps. I'll need to look into #2 later.
-
Mikey,
Thank you for the quick response. I changed MOP(12) to 2 as you suggessted, and according to the new run_gener.csv output file, this did in fact make the generation rates closer to what I was trying to get. However, the most important problem persists. That is, saturations and water volumes in the unsaturated zone don't change through time.
Can anyone help me to figure out why water isn't entering the system?
Mike
-
Can you also check MOP2(3) > 0 with MOP(12)=2?
-
Mike,
Try out the attached input file. I made a few key changes to what you provided:
1. I changed your initial time step size to 1.0e+0 in PARAM.2. It was originally 1.0e10, which forced itough to attempt calculating the whole simulation in a single time step. Since you listed your first time in the TIMES block at 3600. (the end of your injection), it reduced that single step to an hour, rather than 1e+10 seconds. Still, it was unable to resolve the transient behavior of the recharge.
2. I deleted the 3rd record of your PARAM block. Since your initial time step was listed as a positive number, iTOUGH does not know that you're listing individual time steps for it to calculate on the following line.
3. I changed your max time steps (in PARAM.1) from 500 to 9999, which I believe defaults the simulator to look for other stopping conditions (TIMAX or a numerical stability check).
I haven't had a chance to review the outputs too thoroughly, but I definitely see some changes to your output this time.
Let me know how this works.