# Multiple point fault observation association using random forest from analog structural models

Amandine Fratani, Radu S. Stoica, Guillaume Caumon, Jérémie Giraud

<figure><img src="image_1.png" /></figure>

## To cite this version:

Amandine Fratani, Radu S. Stoica, Guillaume Caumon, Jérémie Giraud. Multiple point fault observation association using random forest from analog structural models. 2025. hal-05032543

HAL Id: hal-05032543

https://hal.science/hal-05032543v1

Preprint submitted on 14 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiple point fault observation association using random forest from analog structural models

Amandine Fratani¹,², Radu Stefan Stoica³, Guillaume Caumon¹,⁴, and Jeremie Giraud¹,⁵

¹RING, GeoRessources / ENSG, Université de Lorraine / CNRS, F-54000 Nancy

²ASGA, 2 rue du Doyen M. Roubault, F-54000 Nancy, France

³Université de Lorraine, CNRS, IECL, Inria, F-54000 Nancy, France

⁴Institut Universitaire de Fance (IUF)

⁵Centre for Exploration Targeting (School of Earth Sciences), University of Western Australia, Perth WA 6009, Australia

September 2024

## Abstract

The association of incomplete fault observations is a complex task for which there are several solutions. This problem has already been addressed using a probabilistic approach, where pairwise expert rules are considered. Extending this approach to consider multi-point interactions is an avenue to better represent geological knowledge, but multi-point rules are difficult to define from expert knowledge. To alleviate this limitation, a Random Forest is proposed to infer, from analog three-dimensional models, the potential that a triplet (or a pair) of sparse observations belongs to the same fault. A mathematical formalism is proposed to reconcile the two potentials inside a graph formalism. Two Random Forest learners are trained, one from pair of fault features, the second from triplet of fault features. The features are computed from fault traces extracted from a known 3D geological model (e.g. the length of the fault trace, the throw value, etc.). The learners are tested on fault traces extracted on a contiguous sector of the same model. Triple point potentials are different from existing combination of pair potential, showing that the proposed methodology takes into account multi-point interaction effects.

## 1 Introduction

The presence of underground faults impact subsurface fluid flows which control the formation of strategic mineral resource deposits (e.g. RABEAU et al. (2013)), as well as compartmentalization or drainage in hydrocarbon reservoirs (e.g. BOTTER et al. (2017); LIU et al. (2020); WATTS (1987); ZENG et al. (2021)), CO₂ reservoirs (e.g. NICOL et al. (2017)), and geothermal reservoirs (e.g. FORSTER and EVANS (1991); FORSTER et al. (2003); SILER et al. (2019)). However, identifying fault networks from subsurface data is challenging, in particular for many energy transition applications in which three-dimensional reflection seismic images, using to interpret fault surfaces directly, are not always available. Most often, fault geomodelling relies on “joining the dots” between sparse structural interpretations from geophysical sections, geological sections and boreholes images. On geophysical sections, such as reflection seismic images (FOSSEN, 2016) and geological cross-sections, interpreted faults are represented by identified polygonal lines.

The interpretation of the fault structure requires the association of the fault lines, also named the “Fault correlation problem” by FREEMAN, YIELDING, and BADLEY (1990). In FREEMAN, YIELDING, and BADLEY (1990), the evolution of the fault displacements (of fault traces that an expert wishes to associate) on a strike projection are shown to help the expert for validating the proposed association. Generally, experts tend to propose only one association scenario. However, when different possibilities exist for the interpretation of a dataset, CHAMBERLIN (1890) and TARANTOLA (2006) recommend considering multiple hypotheses (or models) to better assess uncertainties. For this, GODEFROY et al. (2019) propose a probabilistic approach of the fault correlation problem from which several fault networks can be simulated. The determination of probabilities that fault observations belong to the same fault is considered using a graph formalism in which nodes represent fault traces and edges carry the potential that two nodes connected by an edge are associated and belong to
<figure><img src="image_3.png" /><figcaption>Figure 1: Two scenarios where the triplet association potentials is more informative than the average of the pair association potentials. Fault observations are represented using red line, intersecting a black horizon.</figcaption></figure>

the same fault. Several parametric expert rules, obtained from analog reasoning, are used to compute this pair potential, for example by evaluating the consistency between the two nodes and the statistical orientation of the fault datasets in the area. Such rules depend on the geological context and require to be appropriately chosen and parameterized. This graph formalism is interesting to represent the potential that two observations belong to the same fault, but it does not explicitly account for multiple-point potential. Therefore, GODEFROY et al. (2021) compute this multi-point potential from the average of the two-point potentials. However, there is no guarantee that such a computation is predictive. For example, the variations of fault size and displacement along strike may produce large pair association potential, but a low three-point (triplet) association potential if triplet of fault observations does not respect the canonical isolated fault model (WALSH & WATTERSON, 1989) which assumes an elliptical shape with the maximal displacement at the center (Fig 1, scenario 1). Conversely, the triplet potential of three fault observations consistent with this fault model could be higher than just considering the average pair association potential (Fig 1, scenario 2). Fault branchings may also be a significant source of complexity for the determination of multi-point association potential.

Therefore, this paper proposes a model to explicitly include multiple point interactions, starting with three-points, in the determination of the probabilities of associating isolated fault observations. The model is based on the combination of pair potentials and triplet potentials. As defining rules for triplet potential computation is a complex task, supervised machine learning is considered to learn them directly from analog data. This requires the definition and selection of features, corresponding to input variables. Choosing which feature to use is crucial and may vary depending on the context. Several features are proposed in the context of the fault correlation problem in Sect. 2.1 and their importances are evaluated in Sect. 3.4. The learning of individual pairs and triplets potential is performed with a Balanced Random Forest algorithm (CHEN, LIAW, & BREIMAN, 2004) presented in Sect. 2.2. The mathematical model for calculating the association probability is based on a combination of the pair and triplet potentials (Sect. 2.3).
<figure><img src="image_4.png" /><figcaption>Figure 2: Graphic summary of the test case were pair and triplet association potentials estimated from synthetic fault traces using Random Forest. For visualization, only a few pairs and triplets are drawn in the results (d, e and f). f₁ and f₂ correspond to two features used to discriminate the classes.</figcaption></figure>

To evaluate the presented model, synthetic fault traces, called fault observations, are proposed to be generated from three-dimensional model and used in a partly interpreted case (Sect. 3.1), in which the domain is separated in two geographically disjoint subsets: one is used to train the model and the other one to test it (Fig. 2 a and b). Two geological model interpreted by experts are used to generate a simple use case and a more complex one (Sect. 3.2). Two Balanced Random Forests are trained, one with pair features (for example f₁ and f₂ on Fig. 2 c), the other with triplet features, to identify and capture existing associations (Sect. 3.3) from a training database where pairs and triplets of fault observations are identified as either class 1 for associated nodes (belonging to the same fault) and class 0 for disjoint nodes (belonging to different fault surfaces). The performances of the RF trained with different features are presented in Sect. 3.5. For pair association, the resulting potential computed on the test set can be presented directly on the graph (Sect. 3.6 and Fig. 2 d). The triplet potentials, however, cannot be directly visualized on the graph (Fig. 2 e). Therefore, a projected potential on each edge is computed based on the potential of triplets which include this edge (Fig. 2 f). The triplet potential is compared to the average of the three pair potentials as in the work of GODEFROY et al. (2021) to evaluate the interest of considering them (Sect. 3.7). Finally, the proposed model integrating the triplet and pair potentials in a graph is evaluated (Sect. 3.8). The methodology and results are then discussed in Sect. 4.
# 2 Graph formalism for the multipoint fault observation association probabilities determination

## 2.1 Fault observations and association structure

### Fault data interpretation

Faults can be identified either at the surface level directly on outcrops and extrapolated at depth on cross-sections or in the subsurface on borehole images and reflection seismic images. Here, only fault observed on sections (seismic images or cross-sections) are considered. In reflections seismic images (FOSSEN, 2016) and geological cross-sections, interpreted faults are represented by polygonal lines with a range of different characteristics. Parallel sections are used to reduce the number of characteristics attached to each observation. This work considers four characteristics $\mathbf{f}_i$ for each observation $i$, as used classically in fracture and fault characterization studies (e.g., BONNEAU et al., 2013; GOODWIN, AKER, & RØE, 2022): barycenter, length, slope, and maximal vertical displacement (throw).

### Fault observation association structure

Let $\mathcal{V} = \{v_i\}_{1\le i\le I}$ be the set of $I$ incomplete fault observations that are to be associated in order to form a fault network passing through the observed faults. The set of the four characteristics of the fault observations $\mathcal{F} = \{\mathbf{f}_i\}_{1\le i\le I}$ is attached to the set $\mathcal{V}$. The pairwise association status of fault observations is described by the edges $\mathcal{E} = \{e_{ij}\}_{1\le i<j\le I}$ of the graph $G = (\mathcal{V}, \mathcal{E})$. When known or simulated, the association of two fault observations $(v_i, v_j)$ is expressed as a label $L = \{0,1\}$ attached to each edge $e_{ij}$ of $\mathcal{E}$. $l(e_{ij}) = 1$ means the two corresponding fault observations are associated and $l(e_{ij}) = 0$ they are disjoint (i.e., that they belong to two different faults).

Similarly, let's define $\mathcal{T} = \{t_{ijk}\}_{1\le i<j<k\le I}$, the set of triplets of nodes. A triplet $(t_{ijk})$ is considered either associated ($l(t_{ijk}) = 1$) or disjoint ($l(t_{ijk}) = 0$) if at least two observations of $v_i, v_j$ and $v_k$ belongs to two distinct faults.

The potential of associating fault observations are computed on each element of $\mathcal{E}$ and $\mathcal{T}$ from analog data. The analog data are forming a graph $G' = (\mathcal{V}', \mathcal{E}')$ and a set of triplets $\mathcal{T}'$, from which the potentials will be computed as a function of the edge and triplet feature values. For each edge $e_{ij}$ the potential $\alpha_{ij}$ corresponds to the probability that $l(e_{ij}) = 1$ given $\{\mathbf{f}_i, \mathbf{f}_j\}$. Similarly, $\beta_{ijk}$ for triplet $t_{ijk}$ corresponds to the probability that $l(t_{ijk}) = 1$ given $\{\mathbf{f}_i, \mathbf{f}_j, \mathbf{f}_k\}$ the association triplet potential is $\beta_{ijk}$. The methodology to obtain such potentials is described in the next subsection.

### Pair and triplet features

Traditional supervised machine learning calls for choosing appropriate input of features in order to predict the outputs. Choosing the appropriate features is a complex task, especially when spatial data are considered. The vector features of the tuples of observations (pairs or triplets) could be formed by the concatenation of the four characteristics (barycenter, slope, length and maximum throwx) for each fault observation. However, the geographical coordinates used during the training are not the same as those used during prediction. In geostatistics, this issue is circumvented by the intrinsic stationarity hypothesis, which considers statistics to be invariant by translation. Similarly, this work uses a vector between fault observations, hereinafter referred to as distance vector. For pairs, the three-dimensional distance vector ($\overrightarrow{v_iv_j}$ = $(\delta x_{ij}, \delta y_{ij}, \delta z_{ij})$) between the barycenters of the two observations $v_i$ and $v_j$ is considered. Triplets of fault observations $(v_i, v_j$ and $v_k)$ are considered as three pairs, therefore three vectors are computed in the idea to make a triangle: $\overrightarrow{v_iv_j}, \overrightarrow{v_jv_k}, \overrightarrow{v_kv_i}$, with $v_k$ the farthest observation from $v_i$ in xy plans. Additionally, observations are always sorted in lexicographic order (by increasing x, then y coordinate), to ensure that polarity is consistent. The other fault characteristics could be considered as independent from geographical location. However, the evolution of the throw tends to decrease from the fault center to the boundary of the fault (e.g., TORABI and BERG (2011); WALSH and WATTERSON (1989)), so the difference of throw ($\delta t$) between two observations is considered. Similarly, the difference of length ($\delta l$) and of slope ($\delta s$) is used. Each pair $e_{ij}$ or triplet $t_{ijk}$ of fault observations is represented by a feature vector (respectively, $\mathbf{f}_i$ and $\mathbf{f}_{ijk}$) composed of the distance vector, the difference of slope, length and throw of each fault observations. Formally, $\mathbf{f}_{ij} = [\mathbf{f}_i - \mathbf{f}_j]$ and $\mathbf{f}_{ijk} = [\mathbf{f}_i - \mathbf{f}_j, \mathbf{f}_j - \mathbf{f}_k, \mathbf{f}_k - \mathbf{f}_i]$.
## 2.2 Computation of fault observation association potential

The estimation of the association potential for each pair or triplet can be considered as an individual probability estimation task. Considering a pair $e_{ij}$ with a vector of variables $\mathbf{f}_{ij}$ and the label $l(e_{ij})$, the probability estimation task corresponds to the estimation of $P(l(e_{ij}) = 1|\mathbf{f}_{ij})$. The problem is posed in a same way for triplet potential: for a triplet $t_{ijk}$, $P(l(e_{ijk}) = 1|\mathbf{f}_{ijk})$ is the probability that need to be estimated. This problem could, in principle, be treated with logistic regression. However, logistic regression is limited to handle linear relations between input features, whereas non-linearities are expected in principle (as shown in Fig. 1), and will also be shown empirically in Sect. 3.1. Therefore, probability estimation can be done with whatever algorithm consistent for nonparametric regression problem (KRUPPA et al., 2014; MALLEY et al., 2012). This includes for example Random Forest (RF), k-nearest neighbors (k-NN), support vector machines (SVM), and others. These algorithms compute an estimation of probabilities that here are called potentials.

In this work, RF is considered for the estimation of the pair and triplet potentials because this algorithm is known to require few fine-tuning and to perform well across all dimensions (CARUANA, KARAMPATZIAKIS, & YESSENALINA, 2008). RF is a classical algorithm that aggregates the learning and the response of an ensemble of decision trees (BREIMAN, 2001). Decision trees are nonparametric supervised learning methods (BREIMAN et al., 1984) similar to spatial search structures in computational geometry. RF and Decision trees are detailed in Appendix A for the pair association problem.

The number of faults is always limited for a set of fault observations. For example, if only two cross-sections are considered, an observation from one section can be associated with at most one observation from the other section. As a result, there are many more disjoint pairs than associated pairs, leading to an imbalanced dataset. The classic RF algorithm does not integrate the problem of imbalanced training data, which may be source of bias when classes are poorly separated in feature space. In such cases, several possibilities have been exposed in the literature: using cost sensitive trees (KRAWCZYK, WOŹNIAK, & SCHAEFER, 2014) or globally resampling the dataset before applying RF (oversampling the least observed class or undersampling the dominant class) (FERNÁNDEZ et al., 2018). Alternatively, modifications to the algorithm have been proposed, such as using another split criterion (the Hellinger distance (CIESLAK et al., 2012)), or automatically providing a balanced subset to each decision tree (Balanced Random Forest (CHEN, LIAW, & BREIMAN, 2004)), which is the solution retained here after testing the other ones. In a Balanced Random Forest, the pair association potential $\alpha_{ij}$ estimates the probability $P(l(e_{ij}) = 1|\mathbf{f}_{ij})$ of pair association. Similarly, the triplet potential $\beta_{ijk}$ obtained with a Balanced Random Forest trained using a set of triplets and labels estimates $P(l(t_{ijk}) = 1|\mathbf{f}_{ijk})$.

## 2.3 Construction of fault observation association graph

Fault observation associations can be expressed as the existence of edges inside the graph $G$, as explained in Sect. 2.1. Let us define $X$ the random field on the edges corresponding to the possible labeling of the edges ($\mathcal{E}$), and $\Omega$ the set of all possible labelings. The cardinality of $\Omega$ is $L^{\mathcal{E}}$. The most general expression of the distribution of $X$ is given by

$$ P(X(\mathcal{E}) = x) = \varphi_{|\mathcal{E}|}(\mathcal{E}) \left[ \prod_{\{e_{i_1 j_1}\}} \varphi_1(e_{i_1 j_1}) \right] \left[ \prod_{\{e_{i_1 j_1}, e_{i_2 j_2}\}} \varphi_2(e_{i_1 j_1}, e_{i_2 j_2}) \right] \dots \left[ \prod_{\{e_{i_1 j_1}, e_{i_2 j_2}, \dots, e_{i_k j_k}\}} \varphi_k(e_{i_1 j_1}, e_{i_2 j_2}, \dots, e_{i_k j_k}) \right], \quad (1) $$

where the interaction function $\varphi_k$ is applied to all k-tuples $\{e_{i_1 j_1}, e_{i_2 j_2}, \dots, e_{i_k j_k}\}$.

The proposed model hereafter may be considered similar to Cox processes in spatial statistics (MOLLER & WAAGEPETERSEN, 2003). These processes explain clustering of point processes by distributing independent points on the realisation of a random intensity map. Here, the realisation of such a map is given by the computed potentials $D_\alpha = \{\alpha_{ij}\}$ and $D_\beta = \{\beta_{ijk}\}$ obtained for all pairs ($\mathcal{E}$) and triplets ($\mathcal{T}$) from the training data. Conditionally on this map, the labeling model for $X$ based only on functions $\varphi_1$ of first order is independent and of the form

$$ P(X(\mathcal{E}) = x | D_\alpha, D_\beta) = \prod_{\{e_{i_1 j_1}\}} \varphi_1(e_{i_1 j_1} | D_\alpha, D_\beta) \quad (2) $$
is considered in the remainder of this paper.

Several possibilities can be considered to express the function $\varphi(e_{ij}|D_{\alpha}, D_{\beta})$, such as the mean

$$ \varphi(e_{ij}|D_{\alpha}, D_{\beta}) = \frac{\alpha_{ij} + \beta_{ij}}{2}, \qquad (3) $$

the product

$$ \varphi(e_{ij}|D_{\alpha}, D_{\beta}) = \alpha_{ij}\beta_{ij}, \qquad (4) $$

or a weighted mean (with $a, b \in \mathbb{R}^2$),

$$ \varphi(e_{ij}|D_{\alpha}, D_{\beta}) = \frac{a\alpha_{ij} + b\beta_{ij}}{a+b}, \qquad (5) $$

with $\alpha_{ij}$ the RF potential associated to the edge $e_{ij}$ and $\beta_{ij} = \Psi(B_{ij})$ a function of all RF potentials $B_{ij}$ associated to triplets containing the pair $(v_i, v_j)$ as defined in the next paragraph.

One of the problems of considering the triplet contribution is that the set of triplet potentials has a different cardinality from that of $\mathcal{E}$. For visualization and analysis, the choice has been made to transfer the RF potential of the triplets onto the graph edges by considering a weighted mean of the potentials of triplets containing the edge. For each edge $e_{ij}$, a subset $\mathcal{T}_{ij}$ of $\mathcal{T}$ corresponding to the triplets containing the edge $e_{ij}$ is defined. Let $B_{ij} = \{\beta_{ijk}\}_{1<k<1,k \neq i \neq j}$ be the set of triplet potentials attached to $\mathcal{T}_{ij}$. The sum of all the edge potentials of $B_{ij}$ is either equal to zero or positive. If higher than zero, a statistical law $w(t_{ijk})$ is defined by dividing all the potentials of $B_{ij}$ by the sum over $B_{ij}$

$$ w(t_{ijk}) = \frac{\beta_{ijk}}{\sum_{\beta_{ijk} \in B_{ij}} \beta_{ijk}}. \qquad (6) $$

$\beta_{ij}$ is defined as a weighted mean of all the RF potentials of $B_{ij}$ with weights $w(t_{ijk})_{1 \leq k \leq |B_{ij}|}$ as

$$ \beta_{ij} = \Psi(B_{ij}) = \frac{\sum_{k=1}^{|B_{ij}|} w(t_{ijk}) \cdot \beta_{ijk}}{\sum_{k=1}^{|B_{ij}|} w(t_{ijk})}, \qquad (7) $$

with $\beta_{ijk}$ the potential of the triplet $t_{ijk}$. Using Eq. (6), Eq. (7) simplifies to

$$ \beta_{ij} = \frac{\sum_{b=1}^{|B_{ij}|} \beta_{ijb}^2}{\sum_{b=1}^{|B_{ij}|} \beta_{ijb}}. \qquad (8) $$

# 3 Application on a known regional context

## 3.1 Database Generation Method

This section proposes a method for automatic extraction and labeling of synthetic fault observations from existing three-dimensional structural models interpreted and validated by structural geologists (Fig. 3). On each structural model, a set of parallel vertical sections are created perpendicularly to the main axis of the interpreted fault surfaces, and labeled polygonal lines formed by the intersection of the section and the geological boundaries is used to build the training database. These labels are defined by a unique fault identifier. For each intersection line, the fault characteristics described in Sect. 2.1 are computed. Then, all the pairs and triplets are formed from fault lines belonging to different sections. For each pair $e_{ij}$ (respectively, triplet $t_{ijk}$), the reference association label $l(e_{ij})$ (respectively, $l(t_{ijk})$), is set to 1 if the lines i and j (respectively i, j and k) have the same fault identifiers.

The presented methodology is used on part of the model to create the training database, and separately on the other part to test the capacity of the Random Forest to learn from the existing associations and to extend them to the test database. Note that the split between training and test databases is not done by random selection as often in machine learning. Indeed, spatial redundancy exists between close sections, so such a strategy would positively bias the classifier performance. This work prefers testing the extrapolation ability of the classifier to mimic a practical use case.
<figure><img src="image_5.png" /><figcaption>Figure 3: A section extracted from a 3D geological model. On the section, faults are represented in red and horizon lines in blue. On the bottom image, nodes corresponding to each fault line are represented using a red dot and a number.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: The two three-dimensional models. Upper row: the thirteen faults of Model 1 in perspective (a), and top views, showing the training sections (b) and test sections (c). Bottom row: same views for the nine faults of Model 2.</figcaption></figure>

## 3.2 Application Database

Two models are used in the following application, and presented Fig. 4 a and d. For each model, twenty sections are created in the training area (Fig. 4 b and e) and five in the test area (Fig. 4 c and f). Model 1 corresponds to an extensional context with two horizons and thirteen relatively low angle faults, mainly oriented E-W. This
<table><caption>Table 1: Overview of the number of fault observation and elements in the training and test datasets for Model 1 and Model 2.</caption><thead><tr><th>Models</th><th>Dataset</th><th># obs</th><th># pairs</th><th># triplets</th></tr></thead><tbody><tr><td rowspan="2">Model 1</td><td>Training</td><td>190</td><td>17,095</td><td>968,150</td></tr><tr><td>Test</td><td>48</td><td>920</td><td>8,802</td></tr><tr><td rowspan="2">Model 2</td><td>Training</td><td>96</td><td>4,351</td><td>123,737</td></tr><tr><td>Test</td><td>21</td><td>177</td><td>737</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5: Pairwise relationships between some features of Model 1 training (a) and test (b) datasets. Orange lines correspond to the density of associated class, blue one of the disjoint class. The kernel used for the density estimation is Gaussian with a bandwidth selection defined by Scott's Rule (SCOTT, 1992).</figcaption></figure>

model is relatively simple with regard to the association problem as the fault geometry is close to cylindrical. The complexity comes from the lateral termination of some faults. The numbers of observations, pairs and triplets for Model 1 are presented in Table 1. Model 2 is composed of nine faults and two horizons. Faults in this model do not have the same orientation, which is interesting to generate relatively diverse 2D sections in terms of fault slopes and spacing. As for Model 1, the dataset extracted from Model 2 is presented in Table 1.

The relations between the pair features can be explored using density estimation plots of pairwise relationships. Some of the plots are presented in Fig. 5 for Model 1, and in Fig. 6 for Model 2. On Fig. 5, the two classes (disjoint and associated) are superposed on each domain for training and testing sets. The separation between these classes appears as challenging on these two-dimensional views. Moreover, the test set domain is always included in the training domain, expressing a good representativity of the training features domain. The same observations are made on Fig 6 for Model 2. The different cross-plots suggests that pairs of features are insufficient to separate linearly the classes.

## 3.3 Random Forest Parameters

For the estimation of the potentials, the Balanced Random Forest Classifier from the python Imbalanced Learn library (LEMAÎTRE, NOGUEIRA, & ARIDAS, 2017) is chosen. Like the majority of implementations of that algorithm, the important parameters are the number *N* of trees in the forest, the function to measure the quality
<figure><img src="image_8.png" /><figcaption>Figure 6: Pairwise relationships between some features of Model 2 training (a) and test (b) datasets. Orange lines correspond to the density of associated class, blue one of the disjoint class. The kernel used for the density estimation is Gaussian with a bandwidth selection defined by Scott's Rule (SCOTT, 1992).</figcaption></figure>

of a split (the criterion) and the maximum number of features considered during the searching of the best split. Extensive sensitivity studies have been made on the RF parameters to obtain the best choice to combine accuracy and performance. The maximum number of features controls how many features are randomly dropped for each split in a decision tree. Because the number of features is relatively small in this application, all the features are kept during the split. Different number of trees (100, 200, 300, 400 and 500) were tested for the Random Forest. As the results were only slightly better with the larger number of trees, the number of trees was set to 100 to minimize computational cost. Following MINGERS (1989) which claims the equivalence of the Gini impurity and the information gain as split criteria in each decision tree, the Gini impurity was retained.

## 3.4 Analyzing the influence of features

This part evaluates the influence of the features ($\delta x$, $\delta y$, $\delta z$, $\delta l$, $\delta s$ and $\delta t$) on the classification outcome. The construction of RF enables a computation of the participation of each feature to the decision, it is called the feature importance. For each tree, the accumulation of the impurity decreasing at each split is computed. The feature importance is the mean of the accumulations. The feature importance measures the influence of the features on the results.

To complete the evaluation, several RF are run on subset of features, and a metric is used to compute the performance of each RF. For this, the Brier score (BS) (BRIER, 1950) computes the mean squared differences between the true label $l(e_{ij})$ of an edge $e_{ij}$ and the predicted association probability $\alpha_{ij}$ as

$$BS = \frac{1}{|\mathcal{E}|} \sum_{e_{ij} \in \mathcal{E}} (l(e_{ij}) - \alpha_{ij})^2. \qquad (9)$$

Table 2 shows the Brier score computed on pair potential from several scenarios on Model 1 (Fig. 4 c). The chosen features for each scenario are indicated in the corresponding row using the feature importance during the training.

The best result is the scenario 3 considering the differences of distance ($\delta x$, $\delta y$, $\delta z$) and slope ($\delta s$) with a Brier score of 0.019. The computed potentials of these results are presented in Sect. 3.6. However, the second scenario using only the distance shows that the distance features are quite enough for the RF to capture the
<table><caption>Table 2: Brier score on pair probabilities on different scenarios of features from Model 1. Each row represents a scenario and the considered features are the one with the feature importance obtained after the training written in the cell.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="6">Feature importance</th><th rowspan="2">Brier score</th></tr><tr><th>δx</th><th>δy</th><th>δz</th><th>δs</th><th>δl</th><th>δt</th></tr></thead><tbody><tr><td>1</td><td>0.12</td><td>0.88</td><td></td><td></td><td></td><td></td><td>0.130</td></tr><tr><td>2</td><td>0.07</td><td>0.67</td><td>0.26</td><td></td><td></td><td></td><td>0.022</td></tr><tr><td>3</td><td>0.03</td><td>0.59</td><td>0.21</td><td>0.17</td><td></td><td></td><td>0.019</td></tr><tr><td>4</td><td>0.05</td><td>0.64</td><td>0.24</td><td></td><td>0.07</td><td></td><td>0.024</td></tr><tr><td>5</td><td>0.05</td><td>0.65</td><td>0.24</td><td></td><td></td><td>0.06</td><td>0.025</td></tr><tr><td>6</td><td>0.03</td><td>0.57</td><td>0.20</td><td>0.16</td><td>0.04</td><td></td><td>0.021</td></tr><tr><td>7</td><td>0.04</td><td>0.63</td><td>0.24</td><td></td><td>0.04</td><td>0.05</td><td>0.023</td></tr><tr><td>8</td><td>0.02</td><td>0.58</td><td>0.20</td><td>0.16</td><td></td><td>0.04</td><td>0.022</td></tr><tr><td>9</td><td>0.02</td><td>0.58</td><td>0.19</td><td>0.15</td><td>0.03</td><td>0.03</td><td>0.021</td></tr></tbody></table>

interaction between pair of nodes in Model 1. This is supported by the feature importance values. In each scenario the two most important features are $δx$ and $δy$. In this particular model, using the length ($δl$) and the throw ($δt$) does not improve the results, but using the slope $δs$ does. Compared to the feature importances of $δl$ and $δt$, which are always less than 10%, the importance of $δs$ is between 15% and 20% in all scenarios where it is used. This value is slightly below that of the second most important feature. The small feature importance of $δl$ and $δt$ can be explained by the strong coherency of orientation of the faults. Also, the domain size is relatively large with regard to the fault sizes. As a result, many faults are truncated (making the fault line length less relevant), and the throw variability is relatively small. While these feature poorly contribute to the prediction, their use in scenario 9 only entails a relatively small decrease of classifier performance, which shows the good robustness of the RF algorithms when the number of features increases. The contribution of the slope is consistent with the existence of two sets of conjugate faults in Model 1. However, the location of the lines appears as more significant because the two sets of opposite dipping faults are located in distinct subdomains (north and south as visible on Fig. 4, which explains the large importance of $δy$).

<table><caption>Table 3: Brier score on pair probabilities on different scenarios of features from Model 2. Each row represents a scenario and the considered features are the one with the feature importance obtained after the training written in the cell.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="6">Features importances</th><th rowspan="2">Brier score</th></tr><tr><th>δx</th><th>δy</th><th>δz</th><th>δs</th><th>δl</th><th>δt</th></tr></thead><tbody><tr><td>1</td><td>0.85</td><td>0.15</td><td></td><td></td><td></td><td></td><td>0.256</td></tr><tr><td>2</td><td>0.63</td><td>0.11</td><td>0.26</td><td></td><td></td><td></td><td>0.243</td></tr><tr><td>3</td><td>0.58</td><td>0.06</td><td>0.11</td><td>0.25</td><td></td><td></td><td>0.220</td></tr><tr><td>4</td><td>0.58</td><td>0.07</td><td>0.16</td><td></td><td>0.19</td><td></td><td>0.236</td></tr><tr><td>5</td><td>0.57</td><td>0.09</td><td>0.20</td><td></td><td></td><td>0.13</td><td>0.195</td></tr><tr><td>6</td><td>0.55</td><td>0.04</td><td>0.07</td><td>0.23</td><td>0.11</td><td></td><td>0.195</td></tr><tr><td>7</td><td>0.55</td><td>0.04</td><td>0.09</td><td></td><td>0.24</td><td>0.08</td><td>0.188</td></tr><tr><td>8</td><td>0.55</td><td>0.06</td><td>0.14</td><td>0.17</td><td></td><td>0.07</td><td>0.219</td></tr><tr><td>9</td><td>0.54</td><td>0.03</td><td>0.06</td><td>0.23</td><td>0.09</td><td>0.05</td><td>0.191</td></tr></tbody></table>

Table 3 shows the Brier score computed on pair probabilities from several feature usage scenarios on Model 2 (Fig. 4 f). The chosen features for each scenario are indicated in the corresponding row using the feature importance during the training part. Differently from the classification on Model 1, adding features globally decreases the Brier score on Model 2. The best scenario are the 7 and 9, which use most features. In terms of global accuracy, the Brier score is higher than for Model 1, illustrating an overall lower classification performance. This can be explained by the greater complexity of Model 2 in terms of geometry. The presence of curved faults reduces the relative importance of $δx$ and $δy$ as compared to Model 1, and confers a relatively
higher importance to the other features.

## 3.5 Performance

The training and testing of the algorithm were conducted on a laptop computer equipped with 16 GB of RAM and an Intel Core i7 CPU. The training time of the different scenarios for features of Model 1 and Model 2 from Sect. 3.4 are presented in Table 4.

<table><caption>Table 4: Training time in seconds of Balanced Random Forest for each scenario of each models.</caption><thead><tr><th rowspan="3">Scenario</th><th colspan="4">Training times (s)</th></tr><tr><th colspan="2">Model 1</th><th colspan="2">Model 2</th></tr><tr><th>pairs</th><th>triplets</th><th>pairs</th><th>triplets</th></tr></thead><tbody><tr><td>1</td><td>0.64</td><td>17.59</td><td>0.29</td><td>2.93</td></tr><tr><td>2</td><td>0.77</td><td>20.50</td><td>0.37</td><td>3.34</td></tr><tr><td>3</td><td>0.86</td><td>22.40</td><td>0.41</td><td>3.79</td></tr><tr><td>4</td><td>0.92</td><td>23.79</td><td>0.42</td><td>4.13</td></tr><tr><td>5</td><td>0.95</td><td>23.63</td><td>0.40</td><td>4.04</td></tr><tr><td>6</td><td>0.96</td><td>24.43</td><td>0.47</td><td>4.64</td></tr><tr><td>7</td><td>0.98</td><td>24.70</td><td>0.45</td><td>4.56</td></tr><tr><td>8</td><td>1.09</td><td>26.47</td><td>0.47</td><td>4.80</td></tr><tr><td>9</td><td>1.14</td><td>28.49</td><td>0.49</td><td>5.33</td></tr></tbody></table>

Considering the high number of examples in the two training sets (see Table 1), the training time is very low in both cases. For Model 1, the training time for pairs is around one second and less than thirty seconds for triplets. For Model 2, less pairs and triplets are used for the training, leading to smaller computing times.

## 3.6 Potential examples

The presented potential results on Model 1 correspond to the scenario 3 using the distance features ($\delta x$, $\delta y$, $\delta z$) and the slope differences ($\delta s$).

### Pair potential

Pair potentials $\alpha_{ij}$ of the test edges of Model 1 computed using the best scenario of Sect. 3.4 are shown in Fig. 7 a. On the graph, the nodes correspond to the fault observations in map view. Fault observations belonging to the same fault in the model are represented with the same color (Fig. 7 c) as in Model 1 (Fig. 4 a).

South of the domain (Fig. 7), the faults are widely spaced, and the edges with high potential connect exclusively the fault observations belonging to the same fault. To the North, where faults are much closer, there is more deviation from the reference pattern, such as the edges with a potential of nearly 1.0 between green (fault 10) and orange (fault 11) nodes. However, most of the fault observations from the same fault are connected with a high probability.

For two points in a simple case, the Random Forest with only 100 trees using the distance between the barycenters, and the difference of slope is able to capture the existing geometric pattern from the training part and to extend them to an unknown part of the model. A reason for these results is that the faults in the training and test sets have a similar orientation, so the training system shows it ability to capture the geometric features of the data.

### Triplet potential

The graph showing the potentials $\beta_{ij}$ obtained from the computation of the weighted mean (Eq. 7) is presented in Fig. 7 b. In the entire domain, almost all nodes belonging to the same faults are interconnected with a high potential. An exception occurs for the two faults consisting of only two nodes (Fault 2 and Fault 5), which is consistent because triplets including these two nodes should have a low probability of association. The nodes
<figure><img src="image_9.png" /><figcaption>Figure 7: Potential results on the test area of Model 1 in map view showing the fault barycenter as nodes and edges colored by potential value. The nodes are colored according to the fault they belong to.</figcaption></figure>

of faults 1 and 3 are exclusively interconnected when considering high potential values (higher than 0.5). In the rest of the domain, the conclusions are similar to the pair association case: some connections with high potential exist between nodes with different fault identifiers.

## 3.7 Triplet potentials vs pair potentials

The goal of this section is to evaluate whether triplet potentials as evaluated with RF are more predictive than triplet potentials estimated from pair potentials. In GODEFROY et al. (2021), multipoint potentials are approached by considering the mean of pair potentials. For each triplet, the mean of the potentials of the three pairs composing the triplet is computed. It is expected that considering the mean will give worse results than direct triplet potential evaluation, because where a direct estimation gives a zero value, the mean of pairs tend
to give a non zero one. To better take into account the pairs that would have a null potential in the three forming a triplet, the harmonic mean is also considered. The performance of the different methodologies are evaluated using the Brier Score (Eq. 11) in Table 5.

<table><caption>Table 5: Brier score on the different triplet potentials results</caption><thead><tr><th>Model</th><th>RF</th><th>Mean</th><th>Harmonic mean</th></tr></thead><tbody><tr><td>1</td><td>0.003</td><td>0.033</td><td>0.004</td></tr><tr><td>2</td><td>0.070</td><td>0.131</td><td>0.075</td></tr></tbody></table>

For the two models, the best Brier scores are obtained with the RF, but the harmonic mean Brier Score are just slightly above. The results are in line with expectations made in the introduction (Fig. 1).

## 3.8 Graph integrating pairs and triplets

In this section, the model proposed in Sect. 2.3 to integrate the triplet potential to the pair potential on the graph G is evaluated. To assess the benefits of this model, the Brier Score (9) is used to compare the probability to the true class of each edge. The Brier score for the two geological models are shown in Table 6, considering the pair potentials ($\alpha_{ij}$) and the triplet potentials projected on the graph ($\beta_{ij}$) and the conditional probabilities (3-5). The weighted mean (5) is computed with $a = 0.8$ and $b = 0.2$, in order to give more influence to the pair potential than to the triplet potential, to not lose the information when a fault is formed with only two points (e.g., Faults 2 and 5 in Fig. 7 b).

<table><caption>Table 6: Brier score of the different potentials on the edges.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Scenario</th><th colspan="5">Brier Score</th></tr><tr><th>Pair potential</th><th>Triplet potential</th><th>Mean (3)</th><th>Product (4)</th><th>Weighted Mean (5)</th></tr></thead><tbody><tr><td>Model 1</td><td>3</td><td>0.020</td><td>0.018</td><td>0.017</td><td>0.016</td><td>0.018</td></tr><tr><td>Model 2</td><td>9</td><td>0.188</td><td>0.154</td><td>0.159</td><td>0.133</td><td>0.138</td></tr></tbody></table>

For Model 1, the Brier scores considering the triplets are lower than when considering pairs. The difference between the three values is really small (0.01), hence not really significant. However, the smallest values correspond to the triplet potential and the product. Considering only the triplet will lead to the disappearance of the fault composed of only two nodes. For this particular geological model, the product seems the better approach to estimate the fault observation association probabilities on the graph edges using the pair and triplet information.

The Brier scores of the three probabilities integrating the triplet potential to the pair potential for Model 2 are lower than the pair potential. The product yields, again, the smallest Brier score and leads to a significant gain of accuracy as compared to the other options. For the two models, the results are similar and suggest that the product is the best option to jointly represent pair and triplet potentials on the graph edges.

# 4 Discussion

## 4.1 On multi-point association potentials

This paper presented a new methodology to compute the probability of multi-point fault observation association. The RF enables a fast and easy computation of association potentials without a complete fine-tuning. The proposed mathematical model which integrate the triplet potential with the pair potential have better results than just considering pairs and triplets separately.

The presented results show that taking into account the triplets improve the estimation of potential, especially in the complex areas. The proposed formalism for joining the triplet to the pairs on the graph can be extended to enable the creation of a graph reflecting multi-point interactions beyond order three. However, managing the large number of possible configurations with n-point associations would raise computational challenges for the learning algorithm. The graph could be sampled with the algorithm of GODEFROY et al. (2021),
to produce faults networks scenarios by interpoling faults between the associated observations. In addition to
the visual analysis proposed in this work, an interesting avenue to quantitatively validate results could be to
analyze the spatial patterns of the produced realizations.

Another contribution of this paper relates to the methodology to create a training dataset from a reference 3D structural model. This simple approach creates a large number of synthetic fault observations from a 3D geological model separated in two parts for training and testing. It then computes, several thousands of pair and triplet examples for the training of the Random Forest. Importantly, the definition of training and test data as two disjoint sectors is required to test the method, as independent splitting between training and test data would likely lead to overfitting. The process could be enriched by creating sections in other directions (not only parallel) or by using incomplete sections and drillholes to improve representativeness of real underground studies.

## 4.2 On machine learning for inferring the association potentials

Machine learning is considered for replacing the expert rules of GODEFROY et al. (2019) by computing the potential association in the frame of a probability estimation problem. This approach is particularly interesting for the three points where no formal expert rules have yet been proposed. The choice of Balanced Random Forest seems appropriate in the light of the results on pair and triplet potentials, but other classifiers such as boosted trees or neural networks could also be investigated in the future. To evaluate the quality of the classification, the Brier Score used is not perfectly adapted in the case of imbalanced datasets (WALLACE & DAHABREH, 2012) so other metrics version should be envisaged to further analyze the results.

In real world applications, the number of features may vary from one observation to the next. For instance, sometimes on seismic images the throw of an observation cannot be estimated. Moreover, the methodology is tested only on fault observations from synthetic sections. In principle, the graph formalism allows the use of other types of data, such as faults interpreted on borehole images. In this case, the features for the machine learning should be adapted, as borehole images do not provide the fault length. Technically, a mix of origin of fault observations could be integrated in a graph. However, in practice the different features may not correspond to all the types of fault observation, going back to the problem of missing values. Several ways to handle such problem are mentioned in the literature, such as deleting the concerned sample point, replacing missing values by a neutral value such as the mean, median or mode or predicting the missing values (ERIAN, REGALADO, & CONRAD, 2021).

Moreover, as a machine learning algorithm is conditioned on data, it will replicate the concept learned on this data, including the errors. The use case considered in Sect. 3 implies that the algorithm tends to reproduce the concept selected by the expert who interpreted the model. This, therefore may indirectly introduce interpretive bias in the solution. To mitigate this bias, a solution could be to generalize the process to several models at the same time and compute potentials on a complete new set. The set of training model could be obtained from a pool of experts and could also be completed with synthetic models to extend the dataset (WU et al., 2020). In any case, the use of ML raises the question of representativeness of the training set to the problem at hand.

## 4.3 Choice of features and observation uncertainties

For Model 1 (Fig 4 a) which is particularly simple, two of the features have no effect on the potential computation. In a more complex situation such as Model 2 (Fig 4 d), these features improve the discrimination of the label in the RF mainly because the variable orientation and curvature of faults decrease the importance of the ($\delta x, \delta y$) separation. Thus, the other features become more informative. Therefore, considering more features could be a potential enhancement of the methodology to adapt to more geological situations. For example, in the case of synsedimentary faults, the use of only one fault throw value should probably be replaced by a vector of values corresponding to the fault activity (i.e. the movement of an active fault). A side effect of considering more features is to increase the dimensionality of the classification task and possibly to increase the risk of overfitting if the training sample is too small.

In addition, this contribution does not consider the uncertainties on the fault observations. If the level of uncertainty of the existence could be estimated with stochastic approach, such as in the work of TATY MOUKATI et al. (2024) or by considering the fault likelihood of interpreted seismic images (WU et al., 2019), it could be
added as a feature. An additional factor which could be considered is the integration of spatial subdomains in which no faults are present (e.g., intact rock intervals identified on borehole images). Incorporating this information could reduce the number of possible associations. However, its integration into the graph presents challenges, as the nodes represent fault observations rather than the absence thereof.

# 5 Conclusions

This paper presents a new methodology to compute and visualize the potential of associating fault observations. The Balanced Random Forest enables the computation of association potentials from a set of data extracted from three-dimensional geological models deemed representative of the tectonic context. A mathematical model to join the triplet and pair potentials in a graph is proposed and could be extended to simulate fault networks. Experimental results confirm that inferring three-point association potentials produces more accurate results than considering only two-point potentials. Moreover, the mathematical model shows a better results than just considering the pairs or the triplets.

# Acknowledgments

This work was performed through the RING project (https://www.ring-team.org/) at Université de Lorraine. We would like to thank for their support the industrial and academic sponsor of the RING-GOCAD Consortium managed by ASGA. Jeremie Giraud’s research was funded by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101032994. Jeremie Giraud also acknowledges support from the Australian Research Council (grant no. IE240100281 – Integrative geophysics Under cover: getting images of the Unknown Unknowns), the Geological Survey of Western Australia, and the Commonwealth Scientific and Industrial Research Organisation. We acknowledge AspenTech for the SKUA-GOCAD 2019 software and TotalEnergies for the Model 2.

# A Random Forest applied to the pair association problem

Let $\mathcal{C} = \{c_{ij}\} = \{(l(e_{ij}), \mathbf{f}_{ij}), e_{ij} \in \mathcal{E}'\}$ be the set of tuples of training vectors and labels for the pair association problem. To form each decision tree $n$, a subset $\mathcal{C}_n$ is sampled in $\mathcal{C}$ by bootstrap. The space $A_0$ containing $\mathcal{C}_n$ is split in two subspaces $A_1$ and $A_2$, called cells (Fig. 8). The splitting process is repeated for each cell until it contains less than a defined number of observations, or until the ratio of different observations $l$ in the resulting subset falls below a given threshold. Let $H$ be the number of cells at the end of the process. The decision of the split is made by considering an impurity function on each cell. A cell is considered as impure if it contains more than one label inside, and the measure of the impurity is based on the proportions of each label in the cell. For this case, the proportion of label $l = 1$ in any cell $A_{m, 0<m<H}$ of the decision tree is

$$p_m = \frac{1}{|A_m|} \sum_{c_{ij} \in A_m} l(e_{ij}), \qquad (10)$$

and the proportion of label $l = 0$ is $1 - p_m$.

Several impurity functions exist to find the best split at each level of the decision tree for categorical targets, such as the Gini impurity or the entropy. The Gini impurity (BREIMAN et al., 1984) measures the likelihood of a random datum to be misclassified if the given class was randomly chosen according to the class distribution in the cell and is equal to the variance of the label in each cell:

$$f_{gini}(A_m) = p_m(1 - p_m). \qquad (11)$$

The entropy impurity (QUINLAN, 1986) measures the amount of information a feature provides about a class using the Shannon Entropy in each cell (SHANNON, 1948). The expression of this criterion is

$$f_{info\_gain}(A_m) = -p_m \log(p_m). \qquad (12)$$
<figure><img src="image_10.png" /><figcaption>Figure 8: Principle of constructing a decision tree, from left to right.</figcaption></figure>

During the construction of a decision tree, in a Random Forest, the splitting process will only act on a randomly selected subset of the considered features to regularize the solution (BREIMAN, 2001). Finally, the probability $P(l(e_{ij}) = 1|\mathbf{f}_{ij})$ for the element $c_{ij}$ in the tree $r(\mathbf{f}_{ij}, \mathcal{C}_n)$ is $p_m$ if $\mathbf{f}_{ij}$ falls in $A_m$.

The above construction process is performed $N$ times to produce a set (forest) of $N$ decision trees. For a given feature vector $\mathbf{f}_{ij}$, the Random Forest then estimates the probability by averaging over all decision trees the indicator proportion in the final cell of each tree (KRUPPA et al., 2014; MALLEY et al., 2012):

$$r_N(\mathbf{f}_{ij}, \mathcal{C}) = \frac{1}{N} \sum_{n=1}^{N} r(\mathbf{f}_{ij}, \mathcal{C}_n), \qquad (13)$$

with $r(\mathbf{f}_{ij}, \mathcal{C}_n)$ the decision obtained with the decision tree formed with the subset $\mathcal{C}_n \subset \mathcal{C}$.

## References

* BONNEAU F, HENRION V, CAUMON G, RENARD P, & SAUSSE J. (2013, July). A methodology for pseudo-genetic stochastic modeling of discrete fracture networks. *Comput. Geosci.*, 56, 12–22. doi: 10.1016/j.cageo.2013.02.004

* BOTTER C, CARDOZO N, LECOMTE I, ROTEVATN A, & PATON G. (2017, February). The impact of faults and fluid flow on seismic images of a relay ramp over production time. *Pet. Geosci.*, 23(1), 17–28. doi: 10.1144/petgeo2016-027

* BREIMAN L. (2001, October). Random Forests. *Mach. Learn.*, 45(1), 5–32. doi: 10.1023/A:1010933404324

* BREIMAN L, FRIEDMAN J, STONE C, & OLSHEN R. (1984). *Classification and regression trees*. Taylor & Francis.

* BRIER GW. (1950). Verification of forecasts expressed in terms of probability. *Monthly weather review*, 78(1), 1–3.

* CARUANA R, KARAMPATZIAKIS N, & YESSENALINA A. (2008). An empirical evaluation of supervised learning in high dimensions. In *Proceedings of the 25th international conference on Machine learning - ICML '08* (pp. 96–103). Helsinki, Finland: ACM Press. doi: 10.1145/1390156.1390169

* CHAMBERLIN TC. (1890, February). The Method of Multiple Working Hypotheses. *Science*, ns-15(366), 92–96. doi: 10.1126/science.ns-15.366.92

* CHEN C, LIAW A, & BREIMAN L. (2004). Using random forest to learn imbalanced data. *Univ. Calif. Berkeley*, 110(1-12), 24.

* CIESLAK DA, HOENS TR, CHAWLA NV, & KEGELMEYER WP. (2012, January). Hellinger distance decision trees are robust and skew-insensitive. *Data Min. Knowl. Discov.*, 24(1), 136–158. doi: 10.1007/s10618-011-0222-1

* ERIAN KH, REGALADO PH, & CONRAD JM. (2021, June). Missing data handling for machine learning models. *IAES Int. J. Robot. Autom.*, 10(2), 123–132. doi: 10.11591/ijra.v10i2.pp123-132

* FERNÁNDEZ A, GARCÍA S, GALAR M, PRATI RC, KRAWCZYK B, & HERRERA F. (2018). *Learning from Imbalanced Data Sets*. Cham: Springer International Publishing. doi: 10.1007/978-3-319-98074-4

* FORSTER CB, & EVANS JP. (1991). Hydrogeology of thrust faults and crystalline thrust sheets: Results of combined field and modeling studies. *Geophys. Res. Lett.*, 18(5), 979–982. doi: 10.1029/91GL00950
FORSTER CB, JONATHAN, CAINE, SCHULZ S, & NIELSON D. (2003). FAULT ZONE ARCHITECTURE AND FLUID FLOW AN EXAMPLE FROM DIXIE VALLEY, NEVADA.

FOSSEN H. (2016). *Structural geology / Haakon Fossen,...* (Second edition ed.). Cambridge: University Press.

FREEMAN B, YIELDING G, & BADLEY M. (1990, March). Fault correlation during seismic interpretation. *First Break*, 8(3). doi: 10.3997/1365-2397.1990006

GODEFROY G, CAUMON G, LAURENT G, & BONNEAU F. (2019, November). Structural Interpretation of Sparse Fault Data Using Graph Theory and Geological Rules: Fault Data Interpretation. *Math. Geosci.*, 51(8), 1091–1107. doi: 10.1007/s11004-019-09800-0

GODEFROY G, CAUMON G, LAURENT G, & BONNEAU F. (2021, February). Multi-scenario Interpretations From Sparse Fault Evidence Using Graph Theory and Geological Rules. *J. Geophys. Res. Solid Earth*, 126(2). doi: 10.1029/2020JB020022

GOODWIN H, AKER E, & RØE P. (2022, January). Stochastic Modeling of Subseismic Faults Conditioned on Displacement and Orientation Maps. *Math. Geosci.*, 54(1), 207–224. doi: 10.1007/s11004-021-09965-7

KRAWCZYK B, WOŹNIAK M, & SCHAEFER G. (2014, January). Cost-sensitive decision tree ensembles for effective imbalanced classification. *Appl. Soft Comput.*, 14, 554–562. doi: 10.1016/j.asoc.2013.08.014

KRUPPA J, LIU Y, BIAU G, KOHLER M, KÖNIG IR, MALLEY JD, & ZIEGLER A. (2014, July). Probability estimation with machine learning methods for dichotomous and multicategory outcome: Theory. *Biometr. J.*, 56(4), 534–563. doi: 10.1002/bimj.201300068

LEMAÎTRE G, NOGUEIRA F, & ARIDAS CK. (2017). Imbalanced-learn: A python toolbox to tackle the curse of imbalanced datasets in machine learning. *J. Mach. Learn. Res.*, 18(17), 1–5.

LIU G, ZENG L, HAN C, OSTADHASSAN M, LYU W, WANG Q, ... HOU F (2020, January). Natural Fractures in Carbonate Basement Reservoirs of the Jizhong Sub-Basin, Bohai Bay Basin, China: Key Aspects Favoring Oil Production. *Energies*, 13(18), 4635. doi: 10.3390/en13184635

MALLEY JD, KRUPPA J, DASGUPTA A, MALLEY KG, & ZIEGLER A. (2012). Probability machines: Consistent probability estimation using nonparametric learning machines. *Methods Inf. Med.*, 51(1), 74–81. doi: 10.3414/ME00-01-0052

MINGERS J. (1989, March). An empirical comparison of selection measures for decision-tree induction. *Mach. Learn.*, 3(4), 319–342. doi: 10.1007/BF00116837

MOLLER J, & WAAGEPETERSEN RP. (2003). *Statistical Inference and Simulation for Spatial Point Processes*. CRC Press.

NICOL A, SEEBECK H, FIELD B, MCNAMARA D, CHILDS C, CRAIG J, & ROLLAND A. (2017, July). Fault Permeability and CO2 Storage. *Energy Procedia*, 114, 3229–3236. doi: 10.1016/j.egypro.2017.03.1454

QUINLAN JR. (1986, March). Induction of decision trees. *Mach. Learn.*, 1(1), 81–106. doi: 10.1007/BF00116251

RABEAU O, ROYER JJ, JÉBRAK M, & CHEILLETZ A. (2013, October). Log-uniform distribution of gold deposits along major Archean fault zones. *Mineralium Deposita*, 48(7), 817–824. doi: 10.1007/s00126-013-0470-7

SCOTT DW. (1992). *Multivariate Density Estimation: Theory, Practice, and Visualization*. John Wiley & Sons.

SHANNON CE. (1948, July). A mathematical theory of communication. *The Bell System Technical Journal*, 27(3), 379–423. doi: 10.1002/j.1538-7305.1948.tb01338.x

SILER DL, FAULDS JE, HINZ NH, DERING GM, EDWARDS JH, & MAYHEW B. (2019, January). Three-dimensional geologic mapping to assess geothermal potential: Examples from Nevada and Oregon. *Geothermal Energy*, 7(1), 2. doi: 10.1186/s40517-018-0117-0

TARANTOLA A. (2006, August). Popper, Bayes and the inverse problem. *Nature Physics*, 2(8), 492–494. doi: 10.1038/nphys375

TATY MOUKATI F, STOICA RS, BONNEAU F, WU X, & CAUMON G. (2024, July). From Fault Likelihood to Fault Networks: Stochastic Seismic Interpretation Through a Marked Point Process with Interactions. *Math. Geosci.*. doi: 10.1007/s11004-024-10150-9

TORABI A, & BERG SS. (2011, August). Scaling of fault attributes: A review. *Marine and Petroleum Geology*, 28(8), 1444–1460. doi: 10.1016/j.marpetgeo.2011.04.003

WALLACE BC, & DAHABREH IJ. (2012, December). Class Probability Estimates are Unreliable for Imbalanced Data (and How to Fix Them). In *2012 IEEE 12th International Conference on Data Mining* (pp. 695–704). doi: 10.1109/ICDM.2012.115
WALSH JJ, & WATTERSON J. (1989, January). Displacement gradients on fault surfaces. *Journal of Structural Geology*, 11(3), 307–316. doi: 10.1016/0191-8141(89)90070-9

WATTS NL. (1987, November). Theoretical aspects of cap-rock and fault seals for single- and two-phase hydrocarbon columns. *Marine and Petroleum Geology*, 4(4), 274–307. doi: 10.1016/0264-8172(87)90008-0

WU X, GENG Z, SHI Y, PHAM N, FOMEL S, & CAUMON G. (2020, July). Building realistic structure models to train convolutional neural networks for seismic structural interpretation. *GEOPHYSICS*, 85(4), WA27-WA39. doi: 10.1190/geo2019-0375.1

WU X, SHI Y, FOMEL S, LIANG L, ZHANG Q, & YUSIFOV AZ. (2019, November). FaultNet3D: Predicting Fault Probabilities, Strikes, and Dips With a Single Convolutional Neural Network. *IEEE Transactions on Geoscience and Remote Sensing*, 57(11), 9138–9155. doi: 10.1109/TGRS.2019.2925003

ZENG L, LYU W, ZHANG Y, LIU G, & DONG S. (2021). The Effect of Multi-Scale Faults and Fractures on Oil Enrichment and Production in Tight Sandstone Reservoirs: A Case Study in the Southwestern Ordos Basin, China. *Frontiers in Earth Science*, 9.