Building and using fully unstructured (disu) grids

[RyanConway91]

I would like simulate underground mine workings in detail in a regional (aka large model). Seems like a full on disu would be good for this. So far I have stuck to disv.

What are other experience with disu?
What software can be used to generate the grid (I have GWV, which can do octree, but I dont think this will get me quite what I want).
Once you have the grid, intersecting (volumes, point, surfaces) seems non-trivial. How have others done this?

[hjia1005]

I have been using mf6-disu for some complex models. I am happy to share some of my experiences. To create the mesh, you can use:

1. Algomesh (commercial)
2. Gridgen (free and you can create the mesh via Flopy)

For Algomesh, you need to prepare the elevation rasters for each layer and you can set a thickness threshold. Then the extent for each layer can be determined and connections can be set up by Algomesh. You will have a disu file, gsf file, gnode to rnode indexing files, which can be read and processed in Flopy and set up your mf6-disu model in Flopy.

For Gridgen, unfortunately, I have to spend more time on the elevation raster and work out the extent of each layer as polygons beforehand. Then you can use

for layer in nlay:
         g.add_active_domain(f'your_path_to_layer_extent/layer_extent_{layer}', [layer])

By doing this, you get a different number of cells for each layer - that's why we use disu.

After this, you can apply the elevation rasters to your model layers and you need to make your model structure CORRECT before you add any packages - this could involve a lot of iterations of work. The checks may include:

• positive thickness for each cell
• for a specific rnode, top of layer n equals to bottom of layer n-1
• and other checks, until every rnode is correct.

After the above is done you can set up the final mf6-disu file and create a gsf file (this is normally used for pilot point processing because plproc doesn't support mf6-disu yet). You can also create your full mesh as a shapefile and set up your final gnode-rnode file with gnode, rnode, layer, coordinates, top, bottom, thickness for each single cell.

Once you have the mesh in shapefile and the gnode-rnode file (this can be converted to points), you can do all geospatial intersections using geopandas to create other packages (GHB, DRN, RIV, LAK...).

The same way helps to build up a mfusg model, but I personally found that mf6 has a better performance in convergency, especially for a complex model. Also keep in mind that CLN is only available to mfusg, and not for mf6, if you have to use it.

I hope the above gives you some idea to start with.

[RyanConway91]

Thank you this is very helpful. I used the GRIDGEN approach for a test octree DISU and it worked well. Do you use any GUIs to visualize the model? I have so far just exported stuff to vtk and used paraview. I have not used Model Viewer for MODFLOW 6 much (https://www.usgs.gov/software/model-viewer-modflow-6-graphical-user-interface-100-software-release). Is this a good option?

Also, if wanted vertical discretization (for example:)

Do you just use:
g.add_active_domain(f'your_path_to_layer_extent/layer_extent_{layer}', [layer])
To deactivate layers outside of the vertical refinement or is there a better way to directly provide a 3D volume as a refinement feature (instead of a 2D features in algomesh or gridgen) in some other software?

[hjia1005]

I haven't tried Model Viewer yet. I normally export the model to VTK and view it in Paraview. Additionally, I find it useful to export a shapefile containing all layers for viewing it in QGIS.

I have never used vertical discretization yet, if I had to I would define them as new layers and work out their 2D extent beforehand. I still keep the layered concept in all my mf6-disu models to ensure compatibility with other models. In my models, refinement in each layer is consistent, i.e. river cells are refined not only for the top layers but for all layers. This does increase the number of cells in the model, but it's easier for me to create a gnode-rnode table to be used conveniently in subsequent processes.

"Is there a better way to directly provide a 3D volume as a refinement feature?" So far I only use 2D features to set the active domain since my models are literally layered "2D" models. To my current knowledge, I guess you may do this in Paraview by sending the full mesh to Paraview and then intersecting it with a 3D feature. However, after the extraction you will have to work out the connection data to make it recognizable by Modflow - these are actually the dirty parts that FloPy or Algomesh have been dealing with even with the 2D features and they are farrrrrr beyond my capacity. I'll be very interested if you have made any progress or new findings on this topic.

[hjia1005]

Alternatively, you can achieve this without using any additional software or having to redo the connection data. Simply keep the entire model structure. After the 3D intersection in Paraview, set the IDOMAIN of the unwanted cells to 0. This will keep these cells in the model structure, but they will not be included in the simulation. While the model size could be large, the simulation time would remain unchanged.