Divs and disu grid types and ghost nodes

[RyanConway91]

I have been experimenting with disv and disu grids using quad/octree, triangular, and voronoi grid cells. I like quadtree grids best because you can easily make octree disu with gridgen. Do others prefer other grid shapes? If so, what and why?

From the MF6 documentation, quad/octree violate CVFD and ghost nodes should be used. I have not yet done this because I haven't found an easy way to generate the required inputs. Also it sounds like xt3d gets around this? I have tried xt3d and found it makes my models so slow they are basically unusable and also produces very similar results to not using xt3d.

Would running xt3d and observing minimal difference indicate that ghost nodes are not required for a model?

Does anyone have an automated way of generated the gn6 input for quad/oct tree grids?

[javgs-bd]

I tried once Gridgen generated octree grids with different number of nodes per layer, but ended up going back to the quadtree/constant # nodes per layer, because I was getting jagged equipotentials in transitions between refinement zones. I have not had the time of going back and check the cause and have remained intrigued about others experiencies with this.
Looking forward to replies to your questions. Very interesting ones.

[RyanConway91]

Hey @javgs-bd, as I am working through testing different grids, I am noticing the issue with octrees. For examples, here are some simulated concentrations and particle tracks from a disu octree. I am using xt3d and model is converging well:

And here are some mp3du particles:

The concentration clearly shows some strangeness, and for the mp3du run, I had to limit the run time otherwise particles end up bouncing back and forth between cells (my IFACE is correct, the bouncing happens in the middle of the grid where there are no boundary conditions). To me these point something weird in the cell mass balance?

Does anyone know a example octree somewhere or have an idea of what is going on here?

[javgs-bd]

Hi @RyanConway91. I guess my issue was a different one regarding head contours jagged behavior in the transition between grid zones that had different refinement levels. I've never gone into the details of transport in octrees because we went back to use quadtree patch grids.

Although, from what you show, I can only comment on the diamond plume shape. If I recall correctly, this is mentioned in Sorab's BCT manual, and seems to be related with the use of the cross dispersion option (forgive me If I messed up or make up terminology). Of what I'm sure is that there is an example in there with diamond shape plumes.

[RyanConway91]

I am also seeing jagged head contours. I think xt3d is the mf6 equivalent of cross dispersion and is on by default in the .dsp in mf6. Thank you for the discussion on this, I will update if I find the issue.

[hjia1005]

I always use Gridgen to create varying numbers of nodes per layer, but I have never used xt3d and ghost nodes. From my understanding, setting ghost nodes is only necessary when a model fails to converge, correct? My models converge well with reasonable geometry and parameters.

Here are some of the convergence issues I've encountered:

1. Convergence problems sometimes occur with the EVT package when the extinction depth exceeds the thickness of the top cell, if top cells (cells in the top layer, not necessarily Layer 1) are specified as EVT cells. Resetting the depth or assigning more than one cell in a vertical column usually resolves this issue.

2. Convergence issues can also be caused by extremely isolated cells, where a few cells are not connected to other cells or layers in the model. After some trial and error, I mostly avoid these situations by either removing the isolated cells if they are very limited in number and size or trying to improve their connectivity with other cells.

[RyanConway91]

I suppose so. I can see this leading to nasty surprises if you later decide to do a transport or MP3DU run though.

From the stand point of creating a defensible model, seems like the best thing would be to just include gnc6. Im sure even if there were no issue lawyer could use not doing as justification to invalidate the model

[javgs-bd]

The problem is, IMO, that anytime we modellers get near the concept of "valid models", we step into complicated positions.Would then any machine learning or regression models be invalid?
Also, probably, the first and most common GW model violation is that the true principal directions of anisotropy are never aligned with cartesian coordinates. We seldom use (and never know) the full hydraulic conductivity tensor.
I guess in all this is where proper use of model for decision-making based on probabilities and risk come into play.

[hjia1005]

Great discussion and willing to see suggestions from the development team.

[cnicol-gwlogic]

I don't know much about this, but from memory the original mfusg manual had a test case or two where they mapped the errors without gnc. They were tiny from what I recall.

More to the point, I agree with what @javgs-bd said, and further to that - if you're concerned about it, i'd suggest running a comparison case with / without for your particularly worrying case, but not just focusing on "error, whatever that is, but focusing on what implications does it have for your key forecasts of interest - i.e. What the model has ultimately been developed for.

My uneducated guess is that there are far bigger issues at play wrt model utility (and forecast uncertainty (not "error")) than gnc might incur. For the record I've never used gnc or xt3d in my models, although I did play with GNC once and decided not to bother.

[langevin-usgs]

Couple quick thoughts on this:

1. If the CVFD requirement that a connecting line between cell centers must bisect the shared face is violated, then there may be some errors in the resulting solution. The errors are difficult to quantify and notice, because these types of errors do not show up in the mass budget error. These errors also do not seem to cause convergence issues.
2. @cnicol-gwlogic is right in that there is a description of these errors in the MODFLOW-USG manual where there is a comparison between a model with and without ghost nodes. A similar analysis is also presented in the XT3D documentation.
3. With quadtree grids, we recommend smoothing the transition between refinement so that a cell is not connected to more than two cells in any horizontal direction. Gridgen will do this automatically. Errors are larger when grids are not smoothed and refinement changes drastically from one cell to the next.
4. Gridgen also has an option to create the ghost-node information needed to create a GNC package. I don't think this can be activated through flopy, however, but it could be added manually.
5. We generally do not use ghost nodes anymore, because the XT3D approach naturally includes a ghost-node-like correction and it's much easier to use. It's true that XT3D is computationally expensive, but we are working on ways to activate XT3D by connection, so that it is only used where needed.
6. For the most part, if you are using a smoothed quadtree grid, you can probably get away with not using ghost nodes or XT3D, but it would probably be a good idea to do a comparison between a simulation with XT3D or ghost nodes and one without to get a feel for the differences. Echoing what has already been said, for most applications, there are probably larger concerns than these CVFD errors. But keep in mind that these errors can cause wiggles in particle traces and probably additional numerical dispersion in transport models.

[RyanConway91]

@hjia1005, you mentioned you have used Gridgen to create varying numbers of nodes per layer. By this you mean octree/octgrid, correct? If so, do you have issues with "jagged" equipotentials?

For example I am working on 10 layer test model with fine discretization in layer 1 and 10, with the grid smoothed vertically between. I am getting "jagged" heads and concentrations, and MODPATH particles that very erratic. The model converges well.

Have you experienced this before? If so, any tips?

[cnicol-gwlogic]

I think haili was probably talking around uniform mesh but some layers pinching out.

I've heard from others to avoid using vertical refinement, which is what you are talking about, because you often get vertically cyclical fluxes.

[RyanConway91]

By "vertically cyclical flux" you mean a large cell would receive and provide water to smaller cells above or below, creating a vertical loop? That would explain why MODPATH runs indefinitely... Why should this only happen with vertical refinement and not laterally?

So is octree just a grid type that should generally not be used? My model is very simple, just some CHDs and it already seems like its not really working.

If not octree, how would others go about gridding a model of complex surface water (wetlands, streams, ect) over an underground mine? I was hoping to have layer 1 highly refined around the surface water and then have the discretization in deeper layers conform to the mine workings. Seems like I should just use quadtree/triangle and deal with a huge model?

[cnicol-gwlogic]

I'm not an expert on this sorry, but I have always avoided that vertical refinement, despite the big temptation on some painful jobs. I'm not sure why it would happen vertically and not so much laterally, this is just what I have heard second hand.

An alternative might be to use two models in one solution in mf6 - one shallow refined one, another coarser deeper one. Not sure what happens at those exchange boundaries though and whether that is any different from vertical refinement ina single model.

Another alternative may be to reformulate your approach - try to boil the key issue(s) down into more abstract representation(s) - build one or more simple models, each focussed on specific forecast need(s). Maybe then you can avoid the big ugly and painful model! Anyway, just a few thoughts.

Thinking about that a little more, one other thing I have done in your mining / gwsw style example is to save on cells by having lots of layers in and around the mine workings / area of most intense depressurisation (and often stratigraphic mapping detail), and then having coarser layer resolution away from there. You can then connect multiple fine layers to the adjacent coarser layer(s) in the disu. Algomesh has a nice easy method for setting this up ("lateral connection groups"); but you could probably achieve similar with multiple mf6 models in one simulation.

[hjia1005]

Yes I apply the same refinement features and levels across all layers but remove unnecessary cells in some layers. That results in "different number of nodes in different layers"

@cnicol-gwlogic your suggestion above is very insightful, will try to incorporate it in my upcoming models :-)