Adaptive resolution?

[Geomory]

Hello GemPy community,

I am fairly new to GemPy and wanted to create a 3D subsurface model of a project I am working on. A major issue I encountered is that my project area is rather large (e.g. 10 by 10 by 6 kilometers), whereas the structures that I want to depict are rather small. For example, the area is cut by normal faults with geological layers laying near horizontal, subparallel to eachother. The offsets are in the range of 100-400 meters at most. As I understand, by defining the resolution while setting up my model, this discretizes my model space into rectangular boxes with the specified resolution, when I compute the model later on.

Is it somehow possible to refine the model resolution around areas of interest, like faults?
I saw that GemPy offers the functionality of "centered grids" in one of the tutorials about forward gravity.
Could this be used in my case?

Or should I think about cutting down my project area extent in all dimensions and simplifying my dataset?

Thanks in advance and best regards!

[javoha]

Hi @Geomory,
thanks for your question. So first of all you are right that gempy uses dense regular grids by default. This does become a problem in cases like yours, where you want to model small scale features in a large model.
There are three solutions that might solve this:
A) As you said you might need to create smaller models and simplify your dataset. Apart from any grid adaptation the underlying universal-cokriging algorithm will struggle with super complex models, especially if you have many fault interactions. My general suggestion is always to start with the biggest features (main units, dominating faults) and add more things gradually to see where problems arise. It is rarely a good idea to throw a bunch of data for a very complex model at gempy at once. All this however really depens on your specific model, if you could share a skectch or something I might be able to give a little more direction.
B) In the current version of gempy (v3) the feature you are looking for already exists: It's called Octree Grids. This automatically creates grids that are refined around areas of interest (e.g faults, surface boundaries). Right now this feature is kind of in the testing phase and gempy will by default calculate both - an Octree grid and a regular grid. I would have to check if it is possible at the moment to dectivate the regular grid and only calculate on the Octree grid. If not - this is definitely a feature that will be available in the future.
C) While centered grids dont make sense in your case (I think) you can always pass a custom grid. So theoretically you can pass any number of xyz coordinates as a custom grid, deactivate the dense regular grid, and calculate the solution on your custom grid. If you dont want to do the refinement manually, it should also be possible to etract the Octree grid from gempy and use it as a custom grid with a little hacking if B9 doesnt work directly at the moment.

I attached a small image of an example how the Octree grid looks like. I think that would be the optimal solution in your case. Again, if you are able to offer some more details about your model I can also give further advice if you want.

Cheers
Jan

[Geomory]

Hi @javoha,
Thanks for your quick response!
The grid you are showing in your cross section is indeed what I was looking for. I am happy to see that solutions for this exist.
You are right, it is best to start with simple features of my modelling area.

I think I should explain to you a bit more what data I have got and what I am intending to do.

Basically, what I want to do is build a 3D geological model from "a priori" data, that I will describe to you in more detail further below. Once that model is setup, I wanted to introduce uncertainties. For example for the z-values of my geological layers or the x-values on the point data of my faults. This way, I want to generate multiple model realizations. Additionally, my idea was to also calculate the forward gravity of the models and compare them to a real gravity data set that I have for my project area. From what I have seen in the GemPy documentation and the papers that were made by your colleagues, this should be possible to do in GemPy, right?

As for the data and what it looks like:
I have attached a picture of what you can imagine my project area to look like.
Do notice how some faults don't run strictly parallel but are terminated by other faults, this is also the case in my area.

As for my data, I have got:

• Raster data of geological surfaces, which I converted to .csv (X,Y,Z,Formation name).
• Point data of fault locations at several depth interfaces, converted to .csv (X,Y,Z,Fault name).
  In total, I was planning to model 4-5 geological layers and up to 13 faults. But of course, I start simple with the main features and would keep adding complexity.

What I have to do either way is to simplify and reduce the amount of point data that I have available. In total, it is roughly 70,000 points. I learned that GemPy is not meant to handle that much data, and I have to reduce this to be able to compute my model at a reasonable resolution either way. Hence also my initial question about an adaptive "resolution".

The level of detail that I have in my data is somewhat "overkill" for what I want to do, which essentially is stochastic modelling and forward gravity. To understand how I have to simplify my data, I have a few more questions regarding GemPy:

Another point of complexity is that my geological surfaces are not flat. They have a very "low amplitude" topography so to say. You could imagine it like in the picture I attached.

An implication of the topography I have in my layers and the low-offset normal faulting is, that this produces variable offsets along the faults. Is GemPy able to handle something like that? For example, the offset ranges from several 100 meters to several 10 meters along a single fault. Again, this level of detail is way overkill, which is why I thought about simplifying my tectonic blocks by calculating average depths and assuming that they lay horizontally (for now). This way I would have constant offsets along each fault, and I would ensure that it is at least a few hundred meters of offset, so as to not deal with too high of a refinement that I would again run into computational issues.

Do let me know if my concept and my ideas are feasible in GemPy!
At least, from what I have seen so far, I guess it is dependant on the level of complexity I want to keep.

Thanks for your help!

[javoha]

Hi again,
thanks for your detailed description. Let's go through your points:
In general things like these two Graben structures in your sketch can be modelled with gempy. The kind of parallel faults (in depth) that do not crosscut should work well. Gempy also automatically creates a variable offset based on your input data. The deformation (low amplitude folding) will also work, if it is represented in your data.
The main problem I see are truncating faults. Gempy at its current state does not support finite faults (though you can get to an zero offset I guess). And while truncating faults against other faults is possible, this might create issues, especially when creating a complex system and introducing uncertainty.
As suggested before: Start simple and test. Feel free to ask more questions or send a screenshot/piece of code when you encounter issues.
And: Yes - you should definitely downsample your input data.

Cheers.
Jan

[Geomory]

Hi javoha,
Thanks for your inputs and your fast replies!
I wonder, you mentioned that GemPy always computes a regular grid and an octree grid. How does this work when you can only specify either one (resolution or refinement)?
Because as I understand, the refinement level is the number of times the entire 3D space is split up into octants. For example, with an extent of 1000x1000x1000 and a refinement level of 5, the smallest octant will have equal edge lengths of 1000/2^5 ? And then, how does the specified refinement relate to the cell size of the regular grid (and vice versa)?

Also, another question, is there a way to visualize the octree grid together with the computed model like you did in gempy?
I would like to check how the octree grid is refined in my model, but having looked through the documentation, github and tutorials, I have not seen how to do that. Plotting a 2D section works perfectly fine with gempyviewer.plot_2d and the documented arguments.

Also my final question (for now :-)):
As I have a few truncated faults in my model, I wondered if my concept would work at all? I have built a model with artificially generated data, where two normal faults meet at a 90 degree angle (Fault2 is truncated by Fault1, so in the structural frame, Fault1 is younger than Fault2). From my perspective as a beginner in GemPy, it seems to be doing what it should be doing? Except for the fact that Fault2 is also visualized on the other side of Fault1, but I saw that this is a PyVista related thing. I have attached a few screenshots of the artificial model.

Sorry for the walls of text! I very much appreciate your help!

[javoha]

Hi,
so gempy computing both at the moment is kind of a bug not a feature. It has to do with convenience and extraction of the surface meshes that are displayed in 3D visualization. On the one hand a regular grid is the most intuitive grid for new users and also useful for other applications, that's why it is the default. In gempy v3 octree refinement was introduced for faster computation (and it works roughly the way you explained). In gempy version 3 extracting the surface meshes is done by dual contouring which only works on the octree grid. Before marching cubes was used (working on the regular grid).
So right now both are automatically computed (if you dont set a res or refinemen you still get both) but this should be fixed in a future release.

For visualization: You can get the octree grid points:
geo_model.input_transform.apply_inverse(geo_model.solutions.octrees_output[i].grid_corners.values)) (i being the octree level)

But as this is not a mesh you dont get faces. I didn't find a neat way to create and display those - my 2d plot from before was created in a rather hacky way. If you find a nice solution for that I would be happy if you show me. Also, again, the octrees a re kind of beta, I've found some stuff to be not super intuitive.

Regarding you last question:
I have to admit I have not constructed models with too many faults truncating each other. Your example looks good. As before just add one fault at a time and check if the results check out.

Cheers,
Jan

EDIT: Here is an examle of how it can look like: #934