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