Big alternating in/out fluxes along sloped water table boundary
I am using iTOUGH2-ECO2N to investigate upwelling of CO2-laden brine into a shallow aquifer system, but I am having a problem with the initialization of the flow model and would appreciate some help figuring out how to solve it. I am using a 2-D grid, with two layers representing two aquifers that are flat for a while before sloping upward toward the ground surface.
To set the top boundary as a sloped water table (constrained by previously estimated head distributions in the aquifers), I replaced the nodes above the water table elevation with a single water-saturated element that I set to have a large volume, a small nodal distance, and an atmospheric pressure. The left boundary was likewise set to a constant hydrostatic pressure distribution via large-volume, small-distance water-saturated nodes. The other two boundaries were no-flow, as justified by the local geology.
When I run this simulation for a long period of time, I get a pressure distribution that makes sense, and a head distribution (calculated as pressure minus elevation-times-specific-weight-of-water) that also makes sense, given that I'm trying to establish flow from left to right by implementing the sloped water table at the top. These pressure and head distributions are shown in Figures 1 and 2 attached.
However, in order to maintain the pressure conditions that I prescribed, the code calculates significant alternating fluxes into and out of the top boundary, as is shown in Figure 3. While the net water flow rate out of this boundary matches with the total inflow from the left boundary (shown in Figure 4), the alternating fluxes at the top have significant unwanted effects on transport within the system. This is exemplified by the very "jagged" salt plume that forms when I introduce salt and water sources along part of the bottom of the domain, as shown in Figure 5.
The fluxes seem to be related to the head drops that occur at each "stair step" along the sloped top boundary, so I refined the mesh to make these stair steps half their original size, but still got significant alternating fluxes (up to ~25 g/s, as opposed to a maximum of ~30 g/s in the first simulation).
Any help solving this problem would be greatly appreciated!
Thank you in advance,
I'm going to ask a couple of questions that likely won't get to the root of your issue, but I thought it could be a start, so please bear with me. : )
1. As a minor edit, have you considered using dummy elements on the left and upper boundary, instead of using large volumes? This would be an unlikely source of your problems unless you didn't set the volumes of the boundaries large enough.
2. Other than pressure, how are the remaining primary variables set for the left boundary?
3. How are the initial conditions set for the rest of the domain?
I can't promise I'll get to the bottom of your problem, but maybe if we put our heads together we can figure something out. : )
Thanks for your suggestions. Here are my responses:
- Yes, I tried dummy elements and didn't see any change. The 1e50 volume seems to be sufficient to keep the conditions steady.
- On both the top and left boundaries, water saturation is 1 and the salt concentration is set to a constant background value. Temperature is constant throughout the model. I also tried a truly "atmospheric" condition on the top boundary by setting a high gas saturation that would lead to high relative permeability for the gas phase and low relative permeability for the liquid phase, but this caused problems when CO2 started coming into the domain from the top.
- The initial conditions within the domain were much like the boundary conditions; that is, water saturation was set to 1 and salt concentration was set to the constant background value. Pressure inside the domain was given a global initial value of atmospheric, and then I ran TOUGH for a long time to let the pressure distribute according to the boundary conditions, which resulted in the pattern you see in Figure 1.
I appreciate your willingness to help, but I decided to change my approach to avoid the problem I was having. Instead of incorporating the sloped water table boundary, I cut off the domain at the lowest elevation of the boundary, and incorporated constant pressures in large-volume boundary nodes all along the top of the domain. This prevented me from having the "convex" shape of the water table, and from observing what happens above the elevation where I cut off the model, but I'm mostly concerned with the bottom part and the top-right corner anyway.
Please let me know if you think of anything else that's likely to solve the issue in my original model, otherwise I will continue working with the new version.
Thank you very much again for your time,