Exploring Discretization, Grid, and Mesh Types in Groundwater Modeling: Challenges and Experiences

[Groundwater-modeler]

Hi,

I'm not sure if this topic is best suited for the pyemu forum or another platform, but I'll give it a go here since many in the community have likely encountered similar issues.

I have been building models for many years now, but frankly, it can be quite a lonely business when it comes to sharing experiences. Of course, there are vast quantities of reports and materials available online, as well as forums like this one. However, when it comes to cooperative work and practical shared knowledge, I feel that the academic world or specific courses are the best places to exchange experiences and ideas.

Partly because of this, I find the topic I want to discuss a bit difficult to find good sources of information about. I am well aware that the answers I have found in the past (and what I expect here as well) often go along the lines of "it depends on the specific situation," etc.

If you have read this far, please bear with me.

The topic I want to discuss is discretization, grid, and mesh types. First of all, one has to consider which modeling code to use. Modflow and Feflow are codes that come to mind rather quickly. There are probably many others that are more or less known, like HydroGeoSphere and others.

Personally, I prefer Modflow since it is the one I know best and have worked with the most. However, I am intrigued by Feflow's FEM capabilities and its options to create really smooth and well-adapted discretization for all parts of the model, such as fractures, man-made objects, etc.

Modflow's DISU approach seems capable of doing most of what Feflow does well when it comes to adding more freedom to the grid (different grid types, vertical as well as horizontal discretization, and more). I have, for instance, been looking a lot at the octree grid approach to better incorporate dipping fracture zones, mine workings, and tunnel geometries.

From what I understand, though, octree is not supported by pyemu or PEST++. I might have outdated information here. Also, most GUIs do not support octree, although I believe GMS is an exception when it comes to Modflow.

The reason this is interesting is that I feel it would drastically decrease the computational burden of my models and thus make them run smoother and with a more realistic runtime. In my line of work, a typical model would be to simulate the effects (drawdown, water balance, flow patterns) of a tunnel, a mine, or some excavation. These can, of course, often be generalized to better fit in the grid with moderate to high discretization in areas of interest. The vertical discretization is often the big issue since dipping fracture zones often mean a large horizontal area would need to be discretized to realistically capture the geometry across depth.
Most times my models stretch across an area of say 5 by 5 kilometers with a base cell size somewhere around 40 - 160 meters. You quickly realize that in the case of 40 meters and lets say 10 - 15 model layers, one could easily end up with >2 million active cells which starts take a toll on the computational efficiency. Especially when working in a GUI with a heavy burden on working in the model.

What are your experiences, and how do you handle situations like the ones I describe?

How would you, for instance, go about simulating the effects of a tunnel that stretches vertically over a depth of, for example, 80 meters with an equivalent radius of 5 meters? Do you discretize in all directions, adjust conductance alone, or use other tricks?

I would be very interested to hear about your experiences.

Sincerely,

David

[jtwhite79]

This is a tough issue. Ive not used fully 3-D unstructured grids before and I'm pretty intimidated by them TBH. Ive seen several cases where folks got themselves in trouble with 3-D unstructured - yes, it helps keep the computational cost of a single forward run down. But you also need to factor in all the other costs like pre- and post-processing runs, setting up and the modifying the model as needed, being able to quickly and easily inspect/check/verify results in the sense that the model is doing what you want it to, and not just for 1 run, but maybe even an ensemble of runs. For me, its already difficult to quickly and easily detect issues and fix them with structured grid models...so, for me, given my inadequacies with unstructured models, the true cost of a complex 3-D unstructured model is much higher, to the point that the computational savings of the forward run are probably nullified. But there are folks out there who are very skilled with unstructured models so maybe for them, the cost-benefit is worth it...

But if we zoom back, this is really more of an issue of "complexity", as in how complex does the model need to be to provide adequate decision support, where "adequate" is the key word, but one that is very subjective and dependent on many factors. If in your case, you must use 3-D unstructured to represent the key processes and scales, then maybe you need to use that level of complexity...there is a lot of philosophy in this and GMDSI is the place to go for that kind of knowledge.

[p-ortega]

being able to quickly and easily inspect/check/verify results

I had to find a single node to modify a well package in a massive 3D fully unstructured DISU MFUSG. Even with grid intersect from Flopy and all that, to this day, I’m not sure if I got the right node, haha.

[cnicol-gwlogic]

Agree with Jeremy re: trying to keep things as simple as possible to get the answers you need. I try to think about the key receptors, and then the local mode of pressure impact propagation, then work back from there to a minimal setup I need.

With your tunnel example:

• would 2d slice(s) suffice, eg for type case receptors/locations/conditions? And wth that explore the worst and most likely outcomes.
• if you need to represent the entire tunnel length/depth in 3d, can CLNs do the job? Probably far more efficient than discretising a grid around a long dipping fine scale feature like a tunnel. I have done this before, and you can do all sorts of stuff with multiple linear features and interconnections. You could do that in a structured grid too.

Your dipping fracture zone one is hard. I'd try to abstract that down to a simple representation of its potential hydraulic effects as much as possible, remembering that you probably only need to represent such a feature to the level of detail at which you've got data / receptors against which to assess its effects (real observed or potential). But if for example you had mine workings tracking such structures, I think that would be about as hard as this stuff gets. Not sure what the best answer would be there.

[rhugman]

This thread is pretty cool. Y'all should consider turning some of these things into a gmdsi "intro to problem decomposition" or something :)

[Groundwater-modeler]

This is a tough issue. Ive not used fully 3-D unstructured grids before and I'm pretty intimidated by them TBH. Ive seen several cases where folks got themselves in trouble with 3-D unstructured - yes, it helps keep the computational cost of a single forward run down. But you also need to factor in all the other costs like pre- and post-processing runs, setting up and the modifying the model as needed, being able to quickly and easily inspect/check/verify results in the sense that the model is doing what you want it to, and not just for 1 run, but maybe even an ensemble of runs. For me, its already difficult to quickly and easily detect issues and fix them with structured grid models...so, for me, given my inadequacies with unstructured models, the true cost of a complex 3-D unstructured model is much higher, to the point that the computational savings of the forward run are probably nullified. But there are folks out there who are very skilled with unstructured models so maybe for them, the cost-benefit is worth it...

Thanks @jtwhite79 for elaborative and informative response! I agree that there are significant challenges with a fully unstructured approach, especially when one is more accustomed to the traditional i, j, k-structured grids that make matrix navigation inherently simpler. However, I have had positive experiences with DISV-Quadtree, and I don't find the threshold insurmountable when it comes to detecting and addressing issues. This might be because I primarily work with my data through a GUI, which provides a visual representation of the model I'm building, making it easier to identify and rectify problems.

But if we zoom back, this is really more of an issue of "complexity", as in how complex does the model need to be to provide adequate decision support, where "adequate" is the key word, but one that is very subjective and dependent on many factors. If in your case, you must use 3-D unstructured to represent the key processes and scales, then maybe you need to use that level of complexity...there is a lot of philosophy in this and GMDSI is the place to go for that kind of knowledge.

Once again, I agree with almost everything you mentioned. The subjectivity in modeling is indeed a challenge. In an ideal world, multiple models built by different users would produce essentially the same or at least equivalent results. I would like to believe that this is the case, but it may be a naive hope.

Personally, I find it challenging to assess which outcrops or geological materials on a smaller scale are important to represent beforehand, as well as the distance at which they have significant relevance. Since I operate within a commercial context, there is rarely room to build multiple setups of the model and sequentially dismiss various hypotheses.

I can imagine that the geology, as well as the projects, differ significantly between the places we are situated. Here, I would say that the dominant geology consists of topsoil or dry crust on clay overlying till, followed by bedrock. It is often crucial to represent this stratification as the drawdown in the lower confined aquifer is of high interest.

Therefore, if you were to discretize away large areas with bedrock outcrops or where different geological units meet in the horizontal direction, one would lose much of the conceptual understanding of the site I believe.

I guess herein lies the conundrum which often drives me towards higher discretization, often at computational cost...

If I were to reach out to GMDSI, is the GMDSI GitHub the place to go?

Sincerely
David

[Groundwater-modeler]

Agree with Jeremy re: trying to keep things as simple as possible to get the answers you need. I try to think about the key receptors, and then the local mode of pressure impact propagation, then work back from there to a minimal setup I need.

With your tunnel example:

• would 2d slice(s) suffice, eg for type case receptors/locations/conditions? And wth that explore the worst and most likely outcomes.
• if you need to represent the entire tunnel length/depth in 3d, can CLNs do the job? Probably far more efficient than discretising a grid around a long dipping fine scale feature like a tunnel. I have done this before, and you can do all sorts of stuff with multiple linear features and interconnections. You could do that in a structured grid too.

Your dipping fracture zone one is hard. I'd try to abstract that down to a simple representation of its potential hydraulic effects as much as possible, remembering that you probably only need to represent such a feature to the level of detail at which you've got data / receptors against which to assess its effects (real observed or potential). But if for example you had mine workings tracking such structures, I think that would be about as hard as this stuff gets. Not sure what the best answer would be there.

Thank you Chris for your insightful response, I appreciate it.

I agree with your approach of keeping things as simple as possible to get the answers we need. However, finding the right level of simplicity is indeed a challenging task, as I'm sure we are all well aware..

Regarding the use of 2D slices, while they can be effective in some scenarios, often our projects involve large and complex facilities where using 2D sections might not be practical. The end goal is usually to produce maps that show the influence area for the entire complex, and relying solely on 2D slices could result in excessive post-processing work.

CLNs sounds interesting, and I haven't used these before since I have mainly worked with NWT or modflow 6. Typically, I find that regular drain nodes suffice. The challenge, again, is discretization, as conductance can be a difficult parameter to manage. Especially when the task involves estimating the size of the inflow into a facility, and this isn't a known calibration target. There are also other challenges. One being related to converging flows versus radial flows around, for example, a circular installation. If we have a layer that is, say, 50 meters thick (dz) a bit down in the model, it doesn't represent the installation geometry very well. Instead, we get a planar or radial inflow towards these cells rather than a flow converging to a point or a pipe. These types of simplifications are part of the reason why I would prefer to represent installation geometries in my matrix in a more detailed manner. And dont get me started on grouting of tunnels etc.. :)

One could argue that we are in over our heads at some point but when it all comes around, what other options are there? Somehow the best possible decision support predictions must be made and I for one think that you get a lot further with a numerical model than by taking even more conceptual liberties with an analytical counterpart. With numerical modeling, at least parts of the "known" heterogeneity in the system can be represented. With all the flaws that come with it...

Sincerely
David

[neilmanewell]

The timing of this thread is perfect - I was just peer reviewing a Feflow model with 2 million cells, beautiful complexity in the grid to represent the tunnel and stopes (resolution down to centimetres!), and perfect replication of the tunnel/stope progression. The problem was intercepted hydrogeology was represented using one big fat isotropic zone. In my experience, predictive uncertainty reduced by refining the model grid (structural uncertainty) pales in significance to parametric uncertainty. Effort in representing prior hydraulic/recharge uncertainty and the detail in the grid/stresses should be at least equally weighted. All that effort perfectly resolving the tunnels is wasted if the prior parameter distribution is inadequately represented (or in this case, even considered). As others has indicated, there are plenty of resources to aid problem decomposition.

My new ideal maximum runtime for model simulation is less than 5-10 minutes. Work out what you can simplify to achieve these times and start from there. A typical consultant requires these runtimes to enable sufficient parametric freedom, enough agents/realisations to span singular values, and history-match/predict in a reasonable time-frame budget. If that is not possible, then you might have to live with the consequences of incorrectly resolving the posterior (or simply spend $$$$, or use other tools, like DSI/surrogate models etc). Obviously, runtimes like these are harder using finite element solutions, so I favour nimble and numerically stable solutions (3D unstructured MF makes these runtimes achievable).

In summary, all models are crap, and consultants spending countless hours perfectly representing "accuracy" down the scale of metres at the cost of the true likelihood of impacts keeps me up at night.

[jtwhite79]

Agree with that runtime max target!

[Groundwater-modeler]

I couldn't agree more! I usually aim for even lower runtimes, but of course, it is not always achievable. I can very much relate to your description of the Feflow model. Feflow excels in these areas, at least visually. However, I also agree that the overall parameter field is indeed more important than perfect structural representation. For instance, I tend to emphasize representing the bedrock with more layers to accurately depict the declining hydraulic conductivity with depth, coupled with statistical analysis of how this decline behaves.

Regarding your last point, that sounds about right. This is why I am increasingly exploring ensemble-based stochastic approaches to reduce uncertainty, rather than aiming for a single, perfectly calibrated model against some average head measurements over limited time frame in a steady-state scenario. Sometimes modelers spend so much time on calibration that they end up forgetting what the calibration targets actually are and merely accept them as some metric for an absolute truth..

Nevertheless, in cases where it doesn't have to be one or the other, I would favor a relatively accurate representation of the tunnel coupled with a reasonable parameter field over a model that is coarse in both aspects.

[rhugman]

t>1 min makes me itch :)