Correct way to import borehole data

[lzou245]

Hi, I am following this tutorial on transforming borehole data for use within GemPy. As I have had some issues with the subsurface package and my data is relatively simple, I would like to understand how borehole data can be directly transformed into the surface points and orientations used in GemPy.

Here is my code: example.zip

In the input data in the code (an example, the actual format is more complicated), there are a number of straight boreholes containing a coordinate and layer data with the element and the base depth (depth at the bottom of the layer). Based on the tutorial code and what subsurface is doing under the hood, I am creating a StructuralElement for each element, with surface points being the coordinates corresponding to the base depths of each layer with the element present, and no orientation data. Then, I choose an arbitrary point and set its orientation to upwards.

Here is a picture the generated model:

My questions are:

• Is this the correct way to import borehole data into GemPy?
• In the example in the code, element elem3 is not present in the borehole at (0, 0, 100). However, the generated lithological block still extends to (0, 0). I believe I understand why GemPy is doing this, but from a geological perspective I would expect none of elem3 to be present around that borehole because it wasn't measured. Is there a way to implement this behavior with GemPy?
• For real borehole data, there may be multiple layers in the same borehole with the same element that don't belong on a single interface. Is there any way to handle this case?

[javoha]

Hello @lzou245,

thnaks for your question. I will convert this issue to a discussion as I do not think it is a bug or anything similar.

Regarding your three points:

1. While I do not have perosonal experience with the subsurface package this looks correct. Gempy uses points from the bottom of each element as input data and it looks like this is what you achieve here. You can find examples for the minimum data required by GemPy in our new video tutorial series.

2. This is a difficult one. There is no way in gempy to say a layer can absolutely not exist in a certain space. In your example however, you did not specify the order of elements in the structural frame. By default all elements will be in the same strucutral group and therefore feature a conformable behavior (interpolated by a single scalar field) and can by definition not cross-cut. Groupin them into different groups will make them independent and by default younger groups will the erode older groups forming unconfomities. This might lead to elem3 not showing up in the area of that borehole, but it does not guarantee it. I would like to again recommend the video tutorials for a more detailed explanation.
   Here is a screenshot of your model when elem3 is in its own group. It acutally works in this case:

3. Not really. The only solution that would make sense in gempy as far as I can tell would be to rename the "subelements" (e.g sandstaone_a, sandstone_b, or something like this) and treat them as different elements. You can plot them in the sam ecolor to visually acheive the correct effect.

Does this answer your questions?

[lzou245]

Thank you so much for your detailed answer @javoha. Here is my updated example code: example.zip

As per your suggestion to my second question, I am now creating a separate StructuralGroup for each element. For my use case, the element names may be arbitrary so I don't have any prior information on the geological time of each interface. Therefore, I have simply calculated the average z-coordinate for the surface points and ordered the structural groups from highest z to lowest z. In addition, with this new approach, since the surfaces may intersect, they are less useful for me - I will instead focus on the lithological block generated.

I have updated the example data used in my code. The generated block looks like this:

My follow-up questions are:

1. It seems like elem2 is extending to one of the borehole locations without it being present as a layer. I understand that you said even with this approach, it is not guaranteed to disappear around the region. I wanted to confirm again that there is absolutely no way to change this behavior, as it is quite important for me.
2. I noticed that to the right of the red dot near the bottom of the picture, there is a small area seemingly containing elem2. With respect to the generated model, this area being elem2 doesn't make sense to me. Could you explain why this is happening and how I can prevent it?

[lzou245]

Also, I need to do additional processing on the lithological block after it is generated, but I am uncertain about how to extract information about the block from GemPy. I have watched the tutorial video on 3D visualization and export, but it only covers how to export the meshes.

I understand that the block model is stored as a NumPy array in geo_model.solutions.raw_arrays.lith_block. However, consider the below code (works with my updated example above):

lith_block = geo_model.solutions.raw_arrays.lith_block
print(lith_block.shape)
uniques = np.unique(lith_block)
print(uniques)

for group in geo_model.structural_frame.structural_groups:
    print(group.elements[0].id, group.elements[0])

The output is:

(125000,)
[1 2 3 4 5]
2 Element(
        name=elem1,
        color=#9f0052,
        is_active=True
)
1 Element(
        name=elem2,
        color=#015482,
        is_active=True
)
3 Element(
        name=elem3,
        color=#ffbe00,
        is_active=True
)

Firstly, the block clearly needs to reshaped into a 3D array, but what is the order in which the data is specified? My guess is that after doing lith_block.reshape((50, 50, 50)), the (i, j, k)-th entry of the resulting array should correspond to the voxel at (min_x + spacing_x * i, min_y + spacing_y * j, min_z + spacing_z * k) where spacing is extent / resolution. Is this correct?

Secondly, in my example, there are only 3 elements with ids 1, 2, and 3. However, lith_block also contains cells with ids 4 and 5. I am uncertain what this is supposed to represent. Maybe one of them is the basement (?), but why are there two?

Edit: I see that my second question may be related to #1009. Are there any plans for investigating why these additional ids are present, or any possible workarounds?

[javoha]

Hi again,
in regard to your general approach: It is hard to tell if that makes sense for you real model - But the order of the strucutral fram in gempy is very important. Gempy will use this information to mask the result and getting this order worng can really mess with your model.

To your specific questions:

1. I think the only way to brute force this behavir is to set manaul points in these boreholes for the unit that is not supposed to be there above a younger unit. Then the younger element will erode this (if it is in a different structural group). The problem is that this additonal point will influence your interpolation.
2. This is actually a bug. We are working on a fix (as you said in [BUG] Outliers in gempy v3 #1009). I have a rather convoluted cleaning function flying around somewhere but hopefully we will resolve it soonish.
3. The lith_block contains the IDs corresponding to the coordinates in the regular grid.

# Some examples of how to access results
print(geo_model.solutions.raw_arrays.lith_block)
print(geo_model.grid.dense_grid.values)

You can reshape it to the resolution of your grid using

geo_model.solutions.raw_arrays.lith_block.reshape(geo_model.grid.regular_grid.resolution)

The values 4 and 5 in your lith_block both refer to basement (there are only four elements with id 5 which should be a rounding error)

[lzou245]

Thank you @javoha, closing as I do not have any other questions.