Fault relation matrix

[Geomory]

Hi GemPy community,

I am having issues setting up my fault relation matrix as desired. My model consists of 9 structural groups, of which 8 are faults and 1 is my strata. I have set up my geo_model and read the input data for my geo_model before my following code snippet. I am setting up my faults like so:

Create a list of fault file names in the specified order
fault_files = [f"F{i}_simplified.csv" for i in range(1, 9)]

Initialize an empty list to store fault elements
fault_elements = []

Loop through the fault files to read data and create StructuralElement objects

for i, fault_file in enumerate(fault_files, start=1):
    data_fault =pd.read_csv(fault_file)
    
    fault_element = gp.data.StructuralElement(
        name=f"Fault{i}",
        is_active=True,
        color="#2C07FF",
        surface_points=gp.data.SurfacePointsTable.from_arrays(
            x=data_fault["X"].values, 
            y=data_fault["Y"].values, 
            z=data_fault["Z"].values, 
            names=f"Fault{i}"
        ),
        orientations=gp.data.OrientationsTable.initialize_empty()
    )
    
    fault_elements.append(fault_element)

Define a new structural group to contain faults

faults = gp.data.StructuralGroup(
    name="Faults",
    elements=fault_elements,
    structural_relation=gp.data.StackRelationType.FAULT,
    fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_ALL
)

Add the faults group inbetween the sediments group and the default formations.

geo_model.structural_frame.insert_group(1, faults)
print(geo_model.structural_frame)

Add orientations for each fault element

for i in range(1, 9):
    element = geo_model.structural_frame.get_element_by_name(f"Fault{i}")
    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=f"Fault{i}"
    )
print(geo_model.structural_frame)

Map stack to surfaces with all fault groups

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
#        "Sediments": ["Quaternary"]
#    } | {
        "Faults1": ["Fault1"]
    } | {
        "Faults2": ["Fault2"]
    } | {
        "Faults3": ["Fault3"]
    } | {
        "Faults4": ["Fault4"]
    } | {
        "Faults5": ["Fault5"]
    } | {
        "Faults6": ["Fault6"]
    } | {
        "Faults7": ["Fault7"]
    } | {
        "Faults8": ["Fault8"]
    } | {   
        "Strat_Series": ["Tertiary", "Keuper", "Muschelkalk", "Rotliegend"]
    }
)

Set each Fault{i} to fault

for i in range(1, 9):
    gp.set_is_fault(
        frame=geo_model.structural_frame,
        fault_groups=[geo_model.structural_frame.get_group_by_name(f"Faults{i}")]
    )

print(geo_model.structural_frame)

The output is then the following:

print(geo_model.structural_frame)
StructuralFrame(
	structural_groups=[
StructuralGroup(
	name=Faults1,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault1,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults2,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault2,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults3,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault3,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults4,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault4,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults5,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault5,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults6,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault6,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults7,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault7,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults8,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault8,
	color=#527682,
	is_active=True
)
]
),
StructuralGroup(
	name=Strat_Series,
	structural_relation=StackRelationType.ERODE,
	elements=[
Element(
	name=Tertiary,
	color=#728f02,
	is_active=True
),
Element(
	name=Keuper,
	color=#015482,
	is_active=True
),
Element(
	name=Muschelkalk,
	color=#9f0052,
	is_active=True
),
Element(
	name=Rotliegend,
	color=#ffbe00,
	is_active=True
)
]
)
],
	fault_relations=
[[False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False, False]],

However, I want certain faults to be affected by other faults
So, the fault relation matrix I am trying to pass looks like that:

rel_matrix = np.array([
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 1, 0, 0, 1],
  [0, 1, 1, 1, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 1, 1, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 0]])

Now when I run
gp.set_fault_relation(frame=geo_model.structural_frame, rel_matrix=rel_matrix) print(geo_model.structural_frame.fault_relations)

I get the following output:

ValueError: Fault Faults5 cannot affect older fault Faults2

But this can not be correct, since in the structural frame, Faults2 is younger than Fault5?
When I set the corresponding position to "0", it tells me that Faults5 cannot affect older fault Faults3. Same with Faults4.
And when I try to give:

rel_matrix = np.array([
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 1, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],  
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 1, 1, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 0]])

It works, but I get the following fault relationship:

fault_relations=
[[False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False, False, False],
 [False, False, False, False, False, False, False,  True,  True],
 [False, False, False, False, False, False, False, False,  True],
 [False, False, False, False, False, False, False, False, False]],

Which is still not correct. Any pointers to what I am doing wrong?

Thank you!

[javoha]

Hi @Geomory,
thanks for your question. I am sorry that you are having trouble and I think there are maybe some mistakes in the printed output that need fixing. It is also a little hard to find the exact problem just from your snippets but I will try to give some hints.

The general rules for setting this matrix as far as I can tell are:

1. By default faults will cut all older non-fault units (thats what we see in your case, too)
2. Faults can not affect units younger then themselves (thats where I think the Error message is wrong)
3. Faults can not affect themselves

So let's look at your first example:

rel_matrix = np.array([
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 1, 0, 0, 1],
  [0, 1, 1, 1, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 1, 1, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 0]])

The problem here is rule (2) - fault5 ist trying to affect faults 4,3,2 which are all younger - not allowed. You coul however reorder the groups to make fault 5 the youngest and adapt accordingly.

Your second example:

rel_matrix = np.array([
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 1, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],  
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 1, 1, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 1],
  [0, 0, 0, 0, 0, 0, 0, 0, 0]])

Here (and actually also in the example before) as far as I can tell you are trying to set fault7 to cut itself, which is not allowed and therefore overwritten. I am not enitirely sure what is happening with fault4, maybe try setting it to affect everything older than 6 to see if that works.

If you are able to sahre the complete notebook, I might also be able to take a closer look.

Cheers,
Jan

[javoha]

Yes, vonharten(at)cgre.rwth-aachen.de

[javoha]

Hi,
so I at least had time to look into your fault/modeling problem. I was able to reproduce the singular matrix error in your setup.
I checked for duplicate entries (none, which is good) in your input data and applied a nugget effect, but that also did not solve the problem.

Just for reference I computed the model with the default fault relation matrix (bascially faults only affecting the stratigraphic units). Looking at the result leads to a couple of questions for me:

1. It looks like you rock units basically have no offset at fault planes, Do these faults acutally have an effect on the rock units?
2. You have a lot of faults that are cross cutting causing a very complex system of domains. The way the universal co-kriging algorithm handles faults makes the handling of this quite difficult. My guess is that your singular matrix problem resolves from this.

The typical way forward would be: Add faults one by one and always check if the result is as expected. At least that way you might be able to find the core of the problem. You should start with the most dominant features towards the smaller ones.
Personally, I would also advice to leave out faults that do not affect the model in any way for simplification. If you want them in the plot or for other reasons you could add them as meshes to the pyvista plot after model computation.

[Geomory]

Hi,
Thanks for your efforts.

1. Actually they do. But I had to reduce the amount of input points extremely, to be able to compute the model with GemPy at all, and no doubt most of the information on offsets is lost because of this.

2. Yep, this was also my guess that this is the cause of the problem as I encountered it already before I asked for help here. I was able to solve it by tampering with the fault surface points a little, the last time this occured.

I think there must be some sort of bug going on. I already tried adding each fault one-by-one, and even adding only the first three faults and setting the fault matrix accordingly, the fault matrix that is then accepted and output does not match with what I passed.
Using the following structural frame:

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
#        "Sediments": ["Quaternary"]
#    } | {
        "Faults1": ["Fault1"]
    } | {
        "Faults2": ["Fault2"]
    } | {
        "Faults3": ["Fault3"]
    } | {
#        "Faults4": ["Fault4"]
#    } | {
#        "Faults8": ["Fault8"]
#    } | {
#        "Faults5": ["Fault5"]
#    } | {
#        "Faults6": ["Fault6"]
#    } | {
#        "Faults7": ["Fault7"]
#    } | {   
        "Strat_Series": ["Rotliegend"]
    }
)

I try to pass this matrix (fault 3 is cut by fault 2 and all faults cut the strata)

rel_matrix = np.array([
    [0, 0, 0, 1],
    [0, 0, 1, 1],
    [0, 0, 0, 1],
    [0, 0, 0, 1]])
gp.set_fault_relation(frame=geo_model.structural_frame, rel_matrix=rel_matrix)

And the output is:

	fault_relations=
[[False, False, False, False],
 [False, False,  True,  True],
 [False, False, False,  True],
 [False, False, False, False]],

In theory, as I wanted to combine this with the forward gravity method in gempy, all faults should play a role. I wanted to compare the forward model with measured gravity, as a way to validate the structural model (whatever the outcome). But as mentioned earlier, the values are far off the reality and I am not able to reproduce a positive outcome for even a very simple example when using GemPy's forward gravity capabilities.

[javoha]

Hi again,

not sure I get your example completely. The structural frame you show here still has all the faults. I am guessing you just kept 2 faults and the youngest and oldest groups are stratigraphic? The your fault relation matriy, or better the one that gempy generates from your input does make sense, as stratigraphic groups can not offset anything. Maybe I don't get the point though.

Otherwise I feel like even though you say you do have offset, they seem to be very smal compared to the model extent? How did you downsample your points?

[Geomory]

Hey,
No, I did not load the other faults that are commented-out in the map_stack_to_surfaces. The structural frame looks like that:

StructuralFrame(
	structural_groups=[
StructuralGroup(
	name=Faults1,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault1,
	color=#ED0101,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults2,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault2,
	color=#28ED01,
	is_active=True
)
]
),
StructuralGroup(
	name=Faults3,
	structural_relation=StackRelationType.FAULT,
	elements=[
Element(
	name=Fault3,
	color=#F3F700,
	is_active=True
)
]
),
StructuralGroup(
	name=Strat_Series,
	structural_relation=StackRelationType.ERODE,
	elements=[
Element(
	name=Rotliegend,
	color=#00F7F7,
	is_active=True
)
]
)
],
	fault_relations=
[[False, False, False,  True],
 [False, False, False,  True],
 [False, False, False,  True],
 [False, False, False, False]],

And like I said, when I try to make fault 2 cut fault 3, I get

fault_relations=
[[False, False, False, False],
 [False, False,  True,  True],
 [False, False, False,  True],
 [False, False, False, False]],

Yes, the offsets are small compared to the model extent. I hoped I would be able to resolve this sufficiently with the octree grid and as long as I was working with only the faults that do not cut eachother, it worked.

I used vtkQuadricDecimation to downsample the points.