Cross flow between stacked layers
Dears,
Is there is a way in TOUGH3 to output mass flow rate between stacked layers (NOT grid cells) with time?
7 replies
-
I don't know of a direct way to do this in TOUGH3, but you could use COFT on a series of connections between the stacked layers and add up the mass flow rates with a utility program or Python script afterwards. I believe iTOUGH2 does have an option to specify multiple connections and add up the flow rates internally.
-
Just to follow up on Chris's last comment: In iTOUGH2, you can simply specify two material names (or more) to get the total flow across the interface defined by these two materials (i.e., all the connections crossing that interface are identified internally). You could write a relatively straightforward post-processing program that does the same based on the standard TOUGH3 output file, or you directly implement this feature into the TOUGH3 source code.
-
I just released an update to read connection data with toughio (which means that you have to update your previous installation to be able to read connection output data):
pip install -U toughio[full]
I didn't try it myself, but the Python script below should do the trick, assuming:
- Mass flow rates are saved in file OUTPUT_CONNE.csv,
- Elements at the interface are at depth -5.0 m and -15.0 m for layer 1 and layer 2, respectively.
import numpy import toughio # Read MESH file to get element coordinates mesh = toughio.read_input("MESH") # Read connection output data output = toughio.read_output("OUTPUT_CONNE.csv") # Loop over all the connections total_flow = [] for i, (element1, element2) in enumerate(output[0].labels): # Get coordinates of the elements in the connection x1, y1, z1 = mesh["elements"][element1]["center"] x2, y2, z2 = mesh["elements"][element2]["center"] # Check whether the depths are satisfied cond1 = z1 == -5.0 and z2 == -15.0 cond2 = z2 == -5.0 and z1 == -15.0 # Get total flow for current connection for all time steps in OUTPUT_CONNE.csv if cond1 or cond2: total_flow.append([out.data["FLOW"][i] for out in output]) # Convert to numpy array so that rows hold mass flow rate values for all the connections # between the layers and columns represent all the time steps total_flow = numpy.array(total_flow)Note that if your layers are assigned different rock types (e.g. LAY01 and LAY02), you can change the filter condition following
rock1 = mesh["elements"][element1]["material"] rock2 = mesh["elements"][element2]["material"] cond1 = rock1 == "LAY01" and rock2 == "LAY02" cond2 = rock2 == "LAY01" and rock1 == "LAY02"
-
Hi Keurfon,
You're providing excellent services to the TOUGH community - thank you very much!
I have not run your Python script, but it should work for a structured mesh. To make it more generally applicable, I think you would also need to add a test about the direction of a connection across the interface. Recall that the sign of the flow rate is defined by the order of the two elements in a connection, which is arbitrary! For example, if the first element in ALL affected connections is ALWAYS associated with LAY01 and the second with LAY02, everything is fine and you get the correct sum of all flow rates. However, if in some connections the order is reversed, you're adding up flow rates with signs that are just given by a convention rather than the actual flow direction. So check for the "orientation" of the connection with respect to the interface of interest, and flip the sign of the flow rate if needed to make them all consistent. Then also "inform" the user of yoru program, what convention you used, i.e., a positive value means the total flow occurs from LAY02 to LAY01, and negative from LAY01 to LAY02 (which is the convention used in iTOUGH2); this is a bit more difficult to define if the interface is an arbitrary surface, in which case you somehow have to define what is "left" and "right", "front" and "back", "above" and "below", or "in" and "out".
Thanks again!
Stefan
-
Keurfon Luu I meant to respond to this earlier. I've also been working on a Python project to parse TOUGH inputs and outputs similar to (albeit not as well developed as) toughio. In it, I have a boolean-array attribute for each element called "is_n1", which tells if the element is N1 for each of its connections. This would be equivalent to the cond1/cond2 you've put here.