Hecho, Spanish Pyrenees

The Hecho Group turbidites, Spanish Pyrenees

This model is part of a geomodeling benchmarking effort; a full publication describing the benchmark suite is still pending. It is defined on a local, outcrop-scale grid a few meters across, with nine stratigraphic horizons offset by three faults.

import os

import numpy as np

# Aux imports
import pandas as pn

# Importing gempy
import gempy as gp
import gempy_viewer as gpv

Loading surface points from repository:

With pandas we can do it directly from the web and with the right args we can directly tidy the data in gempy style:

data_path = os.path.abspath('../../data/input_data/Hecho')
dfs = []

# First stratigraphic data
for letter in range(1, 10):
    dfs.append(pn.read_csv(
        filepath_or_buffer=data_path + '/H' + str(letter) + '.csv',
        sep=';',
        names=['X', 'Y', 'Z', 'surface', '_'],
        header=0
    ))

# Also faults
for f in range(1, 4):
    fault_df = pn.read_csv(
        filepath_or_buffer=data_path + '/F' + str(f) + 'Line.csv',
        sep=';',
        names=['X', 'Y', 'Z'],
        header=0
    )
    fault_df['surface'] = 'f' + str(f)
    dfs.append(fault_df)

# We put all the surface points together, since that's how gempy likes it:
surface_points = pn.concat(dfs, sort=True)
surface_points.reset_index(inplace=True, drop=False)
surface_points.tail()

	index	X	Y	Z	_	surface
761	6	11.14	-0.17	1.53	NaN	f3
762	7	11.17	0.28	1.68	NaN	f3
763	8	11.16	0.40	1.82	NaN	f3
764	9	11.07	0.03	1.92	NaN	f3
765	10	10.89	-0.43	2.10	NaN	f3

Now we do the same with the orientations:

orientations = pn.read_csv(
    filepath_or_buffer=data_path + '/Dips.csv',
    sep=';',
    names=['X', 'Y', 'Z', 'G_x', 'G_z', '_'],
    header=0
)
# Orientation needs to belong to a surface. This is mainly to categorize to which series belong and to
# use the same color
orientations['surface'] = 0

# We fill the missing direction with a dummy value:
orientations['G_y'] = 0

# Replace -99999.00000 with NaN
orientations.replace(-99999.00000, np.nan, inplace=True)

# Drop irrelevant columns
orientations.drop(columns=['_'], inplace=True)

# Remove rows containing NaN
orientations.dropna(inplace=True)

Data initialization:

Suggested size of the axis-aligned modeling box: Origin: 0 -0.5 0 Maximum: 16 0.5 4.5

Suggested resolution: 0.05m (grid size 321 x 21 x 91)

%%

surface_points_table = gp.data.SurfacePointsTable.from_arrays(
    x=surface_points['X'].values,
    y=surface_points['Y'].values,
    z=surface_points['Z'].values,
    names=surface_points['surface'].values.astype(str)
)

orientations_table = gp.data.OrientationsTable.from_arrays(
    x=orientations['X'].values,
    y=orientations['Y'].values,
    z=orientations['Z'].values,
    G_x=orientations['G_x'].values,
    G_y=orientations['G_y'].values,
    G_z=orientations['G_z'].values,
    names=orientations['surface'].values.astype(str),
    name_id_map=surface_points_table.name_id_map  # ! Make sure that ids and names are shared
)

structural_frame = gp.data.StructuralFrame.from_data_tables(
    surface_points=surface_points_table,
    orientations=orientations_table
)

geo_model = gp.create_geomodel(
    project_name='Hecho',
    extent=[0, 16, -0.5, 0.5, 0, 4.5],
    resolution=[321, 21, 91],
    structural_frame=structural_frame
)

gp.set_section_grid(
    grid=geo_model.grid,
    section_dict={
            'section': ((0., 0.), (16., 0.), (321, 91))
    },
)

Active grids: GridTypes.DENSE|SECTIONS|NONE

	start	stop	resolution	dist
section	(0.0, 0.0)	(16.0, 0.0)	(321, 91)	16.0

We need an orientation per series/fault. Faults do not have an orientation, so the easiest approach is to create one from the surface points available:

f_names = ['f1', 'f2', 'f3']
for fn in f_names:
    element = geo_model.structural_frame.get_element_by_name(fn)
    new_orientations = gp.create_orientations_from_surface_points_coords(
        xyz_coords=element.surface_points.xyz
    )
    gp.add_orientations(
        geo_model=geo_model,
        x=new_orientations.data['X'],
        y=new_orientations.data['Y'],
        z=new_orientations.data['Z'],
        pole_vector=new_orientations.grads,
        elements_names=fn,
        name_id_map=element.surface_points.name_id_map
    )

Now we can see how the data looks so far:

gpv.plot_2d(geo_model)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3e14c7b550>

By default all surfaces belong to one unique series.

geo_model.structural_frame

Structural Groups:	StructuralGroup: Name: default_formation Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: 0 StructuralElement: Name: 0.78 StructuralElement: Name: 1.15 StructuralElement: Name: 1.9 StructuralElement: Name: 2.5 StructuralElement: Name: 3.1 StructuralElement: Name: 3.9 StructuralElement: Name: 4.4 StructuralElement: Name: 5.2 StructuralElement: Name: f1 StructuralElement: Name: f2 StructuralElement: Name: f3
Fault Relations:	default_fo... default_formation
	True False

We need to specify which surface belongs to each series:

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={'Fault1': 'f1', 'Fault2': 'f2', 'Fault3': 'f3'}
)

Structural Groups:	StructuralGroup: Name: Fault1 Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: f1 StructuralGroup: Name: Fault2 Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: f2 StructuralGroup: Name: Fault3 Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: f3 StructuralGroup: Name: default_formation Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: 0 StructuralElement: Name: 0.78 StructuralElement: Name: 1.15 StructuralElement: Name: 1.9 StructuralElement: Name: 2.5 StructuralElement: Name: 3.1 StructuralElement: Name: 3.9 StructuralElement: Name: 4.4 StructuralElement: Name: 5.2
Fault Relations:	Fault1 Fault2 Fault3 default_fo... Fault1 Fault2 Fault3 default_formation
	True False

However, if we want the faults to offset the “Default series”, they will need to be more recent (higher on the pile). We can modify the order by:

Lastly, so far we have not specified which series/faults are actual faults:

gp.set_is_fault(geo_model, ['Fault1', 'Fault2', 'Fault3'])

Structural Groups:	StructuralGroup: Name: Fault1 Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: f1 StructuralGroup: Name: Fault2 Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: f2 StructuralGroup: Name: Fault3 Structural Relation: StackRelationType.FAULT Elements: StructuralElement: Name: f3 StructuralGroup: Name: default_formation Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: 0 StructuralElement: Name: 0.78 StructuralElement: Name: 1.15 StructuralElement: Name: 1.9 StructuralElement: Name: 2.5 StructuralElement: Name: 3.1 StructuralElement: Name: 3.9 StructuralElement: Name: 4.4 StructuralElement: Name: 5.2
Fault Relations:	Fault1 Fault2 Fault3 default_fo... Fault1 Fault2 Fault3 default_formation
	True False

The default range is always the diagonal of the extent. Since in this model data is very close we will need to reduce the range to 5-10% of that value:

geo_model.interpolation_options.kernel_options.range *= 0.2

Explanation of model characteristics and adjustments This model has characteristics that make it difficult to get the right default values: - It is large, and we want high resolution - Some series have a large conditional number (i.e., the model input is not very stable) To address these issues: - Reduce the chunk size during evaluation to trade speed for memory - Reduce the std of the error parameter in octree refinement, which evaluates fewer voxels but may leave some without refinement Enable debugging options to help tune these parameters.

Setting verbose and condition number options for debugging

geo_model.interpolation_options.kernel_options.compute_condition_number = True

gp.compute_model(
    gempy_model=geo_model,
    engine_config=gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.PYTORCH,
        dtype='float64'
    ),
    validate_serialization=True
)

Setting Backend To: AvailableBackends.PYTORCH
GPU enabled. Using device: cuda
GPU device count: 1
Current GPU device: 0
Condition number: 1546449.7334630927.
Chunking done: 21 chunks
Condition number: 1545102.643050272.
Chunking done: 26 chunks
Condition number: 1447785.0114886982.
Chunking done: 21 chunks
Condition number: 40165630.3512226.
Chunking done: 1052 chunks
Chunking done: 9 chunks
Chunking done: 59 chunks
Chunking done: 32 chunks

Solutions: 4 Octree Levels, 12 DualContouringMeshes

gpv.plot_2d(geo_model, cell_number=[10], series_n=3, show_scalar=True)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3e143a8e50>

gpv.plot_2d(geo_model, cell_number=[10], show_data=True)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3e1440cbd0>

gpv.plot_2d(geo_model, section_names=['section'], show_data=True)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3e143a8e50>

sphinx_gallery_thumbnail_number = 3

gpv.plot_3d(geo_model, kwargs_plot_structured_grid={'opacity': .2})

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f3e143efbd0>