Computational modeling and analytical validation of singular geometric effects in fault data using a combinatorial approach

This study analyzes geometrical properties of geological faults using triangulations to model displaced horizons. We investigate two scenarios: one without elevation uncertainties and one with such uncertainties. Through formal mathematical reasoning and computational experiments, we explore how triangular surface data can reveal geometric characteristics of dip-slip faults. Our formal analysis introduces four propositions of increasing generality, demonstrating that, in the absence of elevation errors, duplicate elevation values lead to identical dip directions. For the scenario with elevation uncertainties, we find that the expected dip direction remains consistent with the error-free case. These findings are further supported by computational experiments using a combinatorial algorithm that generates all possible three-element subsets from a given set of points. The results offer insights into predicting fault geometry in data-sparse environments and provide a framework for analyzing directional data in topographic grids with imprecise elevation data.

1 Introduction

Faults influence numerous practical aspects of subsurface geology, including groundwater flow (Bense et al., 2013), hydrocarbon entrapments (Aydin, 2000), geothermal applications (Medici et al., 2023; Fadel et al., 2022), CO₂ storage (Jing et al., 2019), and localized mineralization (Person et al., 2008). In areas with sparse geological data, such as subsurface regions between boreholes, inferring the geometric properties of faults presents significant challenges (Lark et al., 2013; Godefroy et al., 2019). Sparse environments often contain large epistemic uncertainty (uncertainty arising from a lack of knowledge), which can complicate the interpretation of geological structures (Bowden, 2004). While collecting more data can reduce uncertainty (Bond, 2015; Dowd, 2018), practical constraints often make this infeasible.

Recent studies have attempted to reduce epistemic uncertainty in structural geology using triangulations and combinatorial algorithms to analyze fault geometry (Michalak et al., 2021). For example, Michalak et al. (2021) demonstrated that triangles sampled from the walls (see terminology in Fig. 1) of a fault can exhibit counterintuitive behaviors, such as counterintuitive dip directions (towards the upper wall; see Fig. 4b) and identical dip directions from different triangles. This raises intriguing questions about the geometric behavior of triangulated models under sparse data conditions.

This paper builds on this work by providing formal mathematical reasoning that a combinatorial algorithm can reduce epistemic uncertainty in sparse environments. Specifically, here we introduce and analyze the effect of elevation uncertainty on the statistical behavior of the method and present a formal analysis of two scenarios:

1. one with perfect (or rounded) elevation data (Proposition 1, Proposition 2) and

2. one accounting for elevation uncertainties (Proposition 3, Proposition 4).

Then, we explore the statistical implications of our findings using directional data. Following the formal analysis, the work further extends from Michalak et al. (2021) by discussing its relevance to real-world geoscientific datasets. Here, we demonstrate the consequences of these theoretical results in the analysis of 2D and 3D (Fig. 2) directional data derived from topographic grids, which typically consist of points with approximate elevations, commonly observed in bathymetric datasets (Gridded Bathymetry Data, 2024).

2 Background

Uncertainty in geological modeling is a widespread issue affecting various aspects of subsurface analysis. These uncertainties stem from incomplete data, particularly in sparse environments with limited borehole data or surface observations. A key challenge in such cases is accurately relating parts of the study area to lithological units or other geological structures (Wellmann and Regenauer-Lieb, 2012). For example, uncertainty can arise from errors in borehole paths (Pakyuz-Charrier et al., 2018) or in the resulting geological models themselves (Pakyuz-Charrier et al., 2019; Liang et al., 2021). To manage these uncertainties, several methods have been developed, including uncertainty propagation techniques, such as the Markov Chain Monte Carlo (MCMC) method, which estimates uncertainty by feeding model generators with probabilistic input data (de la Varga and Wellmann, 2016; Pakyuz-Charrier et al., 2019).

In sparse geological settings, combinatorial algorithms have emerged as a valuable tool for interpreting fault-related data. One notable example is Godefroy et al. (2019), who developed a method to partition sparse fault evidence into fault clusters using combinatorial techniques. The authors demonstrated that this approach could handle sparse data, but the computational cost increases rapidly as data size grows, governed by Stirling numbers (Allenby and Slomson, 2010).

Orientation measurements, such as dip and dip direction, are critical in subsurface geological modeling and are traditionally collected through fieldwork or outcrop analogs (La Fontaine et al., 2021). These measurements serve as input for co-kriging methods, which combine point data with orientation information to model subsurface structures (de la Varga et al., 2019; Lajaunie et al., 1997; Calcagno et al., 2008). More recently, triangulated data have been widely used in geological modeling to represent 3D surfaces (Merland et al., 2014; Collon et al., 2015; Aydin and Caers, 2017). Triangulations, created by connecting points sampled from geological surfaces, allow the analysis of orientation data by calculating the normal vectors of triangles. This approach is valuable in detecting geological features like faults by clustering these vectors (Michalak et al., 2022).

Despite its utility, triangulation-based analysis faces limitations, particularly when dealing with sparse data environments where the number of triangles available for analysis is reduced. The use of combinatorial algorithms offers a promising alternative by generating all possible triangle configurations from a given dataset (Michalak et al., 2021).

3 Methodology

3.1 Details of the combinatorial approach

We applied the combinatorial algorithm described by Lipski to reduce epistemic uncertainty in determining fault orientations. This algorithm generates all possible three-element subsets (triangles) from a given set of boreholes (Lipski, 2004). This approach systematically creates every possible triangle configuration, enabling a comprehensive geometric analysis. We generated all k-element subsets (k=3) from an n-element set (where n is the total number of borehole locations) to estimate the fault orientation. The algorithm's efficiency allowed us to handle large datasets while ensuring that all potential fault-related triangles were analyzed. The full description of the combinatorial algorithm can be found in Appendix A.

3.2 Singular geometric effects

We begin by recalling and updating the definition of triangles genetically related to faults, originally introduced by Michalak et al. (2021).

Definition 1. A triangle is considered genetically related to a fault (a fault-related triangle) if there exists a pair of its vertices that lie on opposite sides of the fault.

In Michalak et al. (2021), this concept was tied to the notion of non-horizontal triangles. Specifically, a triangle was considered non-horizontal if not all three of its vertices lay on the same side of the fault.

Modification. In this study, due to the introduction of elevation uncertainties, this distinction is no longer sufficient. It is now possible for a triangle to be non-horizontal (in 3D) while all of its vertices still lie on the same side of the fault. Therefore, the definition of a fault-related triangle has been revised accordingly.

Building on the results from Michalak et al. (2021), we extended the analysis of fault-related triangles to account for the observed phenomenon where approximately 8 % of triangles exhibited counterintuitive dip directions (Fig. 4b). We hypothesize that this behavior is controlled by two main factors: the orientation of the edge lying on the horizontal part of the fault (hanging wall or footwall) and the position of the third point relative to this edge (Proposition 1, Fig. 3a). Through formal proofs (see Appendix B), we validated this hypothesis and demonstrated that the dip direction depends on whether the third point lies to the left or right of the fault's edge.

Figure 1 Presentation of the terminology used in the study: (I) the surface of the horizontal footwall, (II) the fault plane, (III) the surface of the horizontal hanging wall, (IV) a horizontal triangle that is genetically unrelated to the fault, (V) a triangle which is genetically related to the fault (a fault-related triangle), (VI) a vertex of the triangle (a point corresponding to a geological interface identified by a borehole). We note that, if elevation uncertainties are introduced, there may be non-horizontal triangles that are genetically unrelated to fault (see Fig. 3d).

3.3 Statistical analysis

Treating the normal vectors as 3D directional data makes it possible to calculate the mean of a group of these 3D vectors. It can be achieved by averaging the Cartesian coordinates of the normal vectors. Then, the resultant vector can be converted to dip direction and dip angle pairs (Allmendinger et al., 2011 – Chap. 2.4). We note that, in this approach, the directional components (X and Y coordinates) are not guaranteed to result in a vector of unit length. Therefore, every vector can contribute differently to the resultant vector (Fig. 2). For example, 3D vectors corresponding to sub-horizontal triangles have smaller values than more inclined triangles at directional components. Therefore, the directional (X and Y) contribution of sub-horizontal triangles to the resultant vector will be relatively small compared to more inclined triangles (Fig. 2).

A different approach would be to conduct a statistical analysis of 2D unit vectors ${\mathit{d}}_{\mathrm{1}},\mathrm{\dots },{\mathit{d}}_n$ corresponding to the initially collected 3D unit normal vectors of triangles $t_{\mathrm{1}},\mathrm{\dots },t_n$, where n denotes the number of observations. The mean direction $\stackrel{\mathrm{‾}}{\mathit{\theta }}$ of 2D unit vectors ${\mathit{d}}_{\mathrm{1}},\mathrm{\dots },{\mathit{d}}_n$ and their corresponding angles ${\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_n$ is defined as the direction of the resultant vector ${\mathit{d}}_{\mathrm{1}}+\mathrm{\dots }+{\mathit{d}}_n$ (Mardia and Jupp, 2008). Firstly, the Cartesian coordinates of the center of the mass (Mardia and Jupp, 2008) are calculated as follows:

$\stackrel{\mathrm{‾}}{C}=\frac{\mathrm{1}}{n}{\sum }_{j=\mathrm{1}}^n\cos {\mathit{\theta }}_j$, $\stackrel{\mathrm{‾}}{S}=\frac{\mathrm{1}}{n}{\sum }_{j=\mathrm{1}}^n\sin {\mathit{\theta }}_j$. We note that, in our case, the X and Y axes are aligned with the north and east directions, respectively. Therefore, the $\stackrel{\mathrm{‾}}{C}$ and $\stackrel{\mathrm{‾}}{S}$ values correspond to the north and east directions, respectively (a different convention is adopted in the textbook of Mardia and Jupp, 2008). To calculate the mean direction $\stackrel{\mathrm{‾}}{\mathit{\theta }}$, we use the following formula (modified from Fisher, 1993):

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{\mathit{\theta }}=\left\{\begin{array}{ll}\text{atan}\left(\frac{\stackrel{\mathrm{‾}}{S}}{\stackrel{\mathrm{‾}}{C}}\right),& \stackrel{\mathrm{‾}}{S}>\mathrm{0},\stackrel{\mathrm{‾}}{C}>\mathrm{0}\\ \text{atan}\left(\frac{\stackrel{\mathrm{‾}}{S}}{\stackrel{\mathrm{‾}}{C}}\right)+\mathit{\pi },& \stackrel{\mathrm{‾}}{C}<\mathrm{0}\\ \text{atan}\left(\frac{\stackrel{\mathrm{‾}}{S}}{\stackrel{\mathrm{‾}}{C}}\right)+\mathrm{2}\mathit{\pi },& \stackrel{\mathrm{‾}}{S}<\mathrm{0},\stackrel{\mathrm{‾}}{C}>\mathrm{0}.\end{array}\right.\end{array}$$

The resultant length is the length of the resultant vector sum ${\mathit{d}}_{\mathrm{1}}+\mathrm{\dots }+{\mathit{d}}_n$, and the mean resultant length is defined as the length of the center of the mass vector $\stackrel{\mathrm{‾}}{R}=\sqrt{{\stackrel{\mathrm{‾}}{C}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{S}}^{\mathrm{2}}}$. We also calculated the median direction using the circular package (Lund et al., 2017), and it is any angle ϕ such that (Mardia and Jupp, 2008)

1. half of the data points lie in the arc $\left[\mathit{\varphi },\mathit{\varphi }+\mathit{\pi }\right]$

2. the majority of the data points are nearer to ϕ than to ϕ+π.

The sample circular standard deviation is defined as $\sqrt{-\mathrm{2}\log \left(\mathrm{1}-V\right)}=\sqrt{-\mathrm{2}\log \stackrel{\mathrm{‾}}{R}}$, where $V=\mathrm{1}-\stackrel{\mathrm{‾}}{R}$ denotes the sample circular variance (Mardia and Jupp, 2008). Using $\mathrm{1}-\cos \left(\mathit{\theta }-\mathit{\xi }\right)$ as a measure of the distance between angles θ and ξ, it can be shown that V can be used as a measure of dispersion around the mean dip direction and that it is equal to $V=D\left(\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)=\frac{\mathrm{1}}{n}{\sum }_{j=\mathrm{1}}^n\mathit{\{}\mathrm{1}-\cos \left({\mathit{\theta }}_i-\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)\mathit{\}}$. We also calculate the sample circular dispersion (Fisher, 1993) defined as $\stackrel{\mathrm{‾}}{\mathit{\delta }}=\frac{\mathrm{1}-m_{\mathrm{2}}}{\left(\mathrm{2}\stackrel{\mathrm{‾}}{R^{\mathrm{2}}}\right)}$, where m₂ denotes the second central trigonometric moment and is equal to $m_{\mathrm{2}}=\frac{\mathrm{1}}{n}{\sum }_{j=\mathrm{1}}^n\cos \mathrm{2}\left({\mathit{\theta }}_i-\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)$. We also use the circular standard error defined as the square root of the sample circular dispersion divided by the number of samples $\stackrel{\mathrm{‾}}{{\mathit{\sigma }}^{\mathrm{2}}}=\frac{\stackrel{\mathrm{‾}}{\mathit{\delta }}}{n}$. Using the above quantities, non-parametric methods can be used to estimate 100(1−α) % confidence intervals for $\stackrel{\mathrm{‾}}{\mathit{\theta }}$:

($\stackrel{\mathrm{‾}}{\mathit{\theta }}-\arcsin \left(z_{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\alpha }}\stackrel{\mathrm{‾}}{\mathit{\sigma }}\right)$, $\stackrel{\mathrm{‾}}{\mathit{\theta }}+\arcsin \left(z_{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\alpha }}\stackrel{\mathrm{‾}}{\mathit{\sigma }}\right))$, where z_½α is the upper $\mathrm{100}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\alpha }\right)\mathit{\%}$ point of the N(0,1) distribution.

Figure 2 Illustration of the difference between calculating the mean direction using x and y components of the 3D and 2D unit vectors. In the first approach, the vectors impact the resulting dip directions differently. Triangles with greater dip angles have a more significant impact, and more horizontally oriented triangles have minor x and y components, making their contribution less important. In the 2D analysis, all vectors are equally important because they all have unit length.

4 Propositions

4.1 A guide to propositions

This section summarizes the main theoretical findings from the four key propositions, focusing on their implications for fault-related triangles (the workflow of the theoretical part is presented in Fig. 3).

Proposition 1 shows that the dip direction of fault-related triangles is controlled by the orientation of the fixed edge on the footwall and the position of the third point on the hanging wall given constant elevation difference between the hanging wall and the footwall. The dip direction can only be one of two values, differing by 180°, and is orthogonal to the fixed edge. The third point's position (left or right of the edge) determines which direction is observed (Fig. 4).

Proposition 2 extends this by proving that adding horizontal triangles from the footwall does not alter the mean dip direction of the fault-related triangles. Since the normal vectors of footwall triangles are $\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]$, they do not influence the dip direction.

Proposition 3 considers a fixed edge on the footwall and a third point on the hanging wall. It investigates the effect of elevation errors and concludes that the expected dip direction remains unchanged, even when elevation uncertainties are introduced. This demonstrates that fault analysis remains robust despite moderate elevation inaccuracies.

Proposition 4 generalizes Proposition 2 for cases with elevation uncertainties. Adding triangles with uncertain elevations does not affect the mean dip direction because their X and Y components average to zero.

Together, these propositions show that fault-related triangles reliably indicate dip direction, even in the presence of elevation uncertainties.

Figure 3 Workflow in the theoretical part of the study. Panel (a) relates to Proposition 1. In this proposition, an edge e is fixed on the surface of the footwall. Then, the dip directions of triangles depend on whether the third point sampled from the hanging wall lies either to the left or to the right of e. The points on the surfaces of the hanging wall and footwall have constant elevation. Panel (b) relates to Proposition 2, and it shows a more general scenario than presented in panel (a). Here, for a fixed edge e, dip directions of triangles depend on whether the third point lies either to the left or right of e. We do not require that the points lie on the hanging wall because the $\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]$ vectors from the footwall will not affect the average dip direction of triangles. (c) Here, we add uncertainty to the elevation of points (Proposition 3). Then, for a fixed edge e, we investigate the average dip direction of triangles with the third point from the hanging wall. (d) The last scenario is the most general scenario (Proposition 4), in which points have uncertain elevations and are no longer required to be sampled from the surface of the hanging wall. In other words, for a fixed edge e, the dip directions of triangles depend on whether the third point lies either to the left or to the right of e. We no longer require that points lie only on the hanging wall. This is because the expected values of X and Y coordinates associated with triangles on the footwall are zero.

4.2 Propositions

In the following analysis, we assume that the considered triangles are non-vertical and non-horizontal. The reason is that horizontal triangles give no information about the dip direction, and, in the case of vertical triangles, there is no possibility of deciding which of the two directions corresponds to the direction of the dip.

Proposition 1.

Let T be a set of non-vertical and non-horizontal triangles genetically related to the fault, and let e:={p₁p₂} be the fixed edge lying on the footwall. When the difference between the hanging wall and the footwall is constant, then the following facts hold:

• (A)

  There are only two possible dip directions for triangles in T, d₁ and d₂, which have a dip direction difference of 180°.

• (B)

  The two different dip directions d₁ and d₂ are orthogonal to e.

• (C)

  Moreover, the specific value of the direction d₁ or d₂ depends on whether the third point p₃ located at the hanging wall lies to the left or right of the line containing e.

Proof. See Appendix B and the corresponding illustration (Fig. 4).

Figure 4 Illustration of Proposition 1. There are three triangles that share the same edge on the footwall. Two of the triangles (yellow and orange) have the same dip direction (302°). In other words, the projections of the yellow and orange vectors onto the horizontal plane are parallel. We note that l(e) denotes the line containing e. Then, the remaining points forming the yellow and orange triangles lie to the left to l(e). The third triangle (cyan) has the opposite direction (122°), and the third point lies to the right of l(e). In panel (c), we present the three orientation measurements from panels (a) and (b) on the spherical projection. The spherical projection was performed using the Stereonet software (Allmendinger et al., 2011; Cardozo and Allmendinger, 2013).

Observation 1. For a fixed horizontal edge on the footwall and a vertical fault dipping to the west (azimuth 270°), the probability that a triangle will dip precisely to the east (azimuth 90°) is zero.

Proof. For a triangle to dip exactly to the east, one should have an edge aligned with the N–S direction. However, the third point must be to the right of such an edge. However, the hanging wall is to the left of this edge, so no appropriate points can be sampled from the surface of the hanging wall or the footwall.

Proposition 2. Extending Proposition 1, we demonstrate that adding horizontal triangles from the footwall does not alter the mean dip direction of fault-related triangles.

Proof. The set of triangles on the right side of the edge can be divided into those with the third point on the hanging wall (indices from 1 to k) and those on the footwall (vertices from k+1 to n).

$$\begin{array}{}\text{(2)}& \begin{array}{l}{\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}\right]\\ {\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}}\right]\\ \mathrm{⋮}\\ {\mathit{v}}_k=\left[x_k,y_k,z_k\right]\\ {\mathit{v}}_{k+\mathbf{1}}=\left[x_{k+\mathrm{1}},y_{k+\mathrm{1}},z_{k+\mathrm{1}}\right]=\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]\\ \mathrm{⋮}\\ {\mathit{v}}_n=\left[x_n,y_n,z_n\right]=\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]\end{array}\end{array}$$

However, we know that all vectors from the footwall starting from v_k+1 to v_n have the form $\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]$ and that they contribute nothing to the dip direction.

Next, we can try to adapt and rewrite Proposition 1 (the third point on the surface of the hanging wall) for uncertain elevation data. It means that we add normally distributed errors (with an expected value equal to zero) to the elevations. Using this assumption, Proposition 1 can be rewritten in the following way.

Proposition 3.

When introducing elevation uncertainties, we find that the expected dip direction remains consistent with the error-free case. The third point is required to be from the surface of the hanging wall.

Proof. See Appendix C.

Proposition 4.

This proposition generalizes Proposition 3, showing that adding footwall triangles with uncertain elevations still does not affect the mean dip direction. The X and Y components of the footwall triangles' normal vectors average to zero, further confirming the robustness of the method under elevation uncertainties.

Proof. As in the deterministic example (Proposition 2), the set of triangles on the right side of the edge can be divided into those with the third point on the hanging wall (indices from 1 to k) and those with the third point on the footwall (vertices from k+1 to n). To complete the proof, we need to show that the expected values of the X and Y coordinates of the normal vectors representing triangles with third points lying on the surface of the footwall (from v_k+1 to v_n) are zero. We present the desired proofs in Appendix D.

$$\begin{array}{}\text{(3)}& \begin{array}{l}{\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}\right]\\ {\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}}\right]\\ \mathrm{⋮}\\ {\mathit{v}}_k=\left[x_k,y_k,z_k\right]\\ {\mathit{v}}_{k+\mathbf{1}}=\left[x_{k+\mathrm{1}},y_{k+\mathrm{1}},z_{k+\mathrm{1}}\right]\\ \mathrm{⋮}\\ {\mathit{v}}_n=\left[x_n,y_n,z_n\right]\end{array}\end{array}$$

This implies that the sums of the X and Y coordinates from v_k+1 to v_n are zero: ${\sum }_{k+\mathrm{1}}^n{x}_i={\sum }_{k+\mathrm{1}}^n{y}_i=\mathrm{0}$. Therefore, as in the deterministic case, adding these sums will not change the first two coordinates of the expected normal vector. Appendix D provides further discussion on when the mean vector from v_k+1 to v_n is parallel to $\left[\mathrm{0},\mathrm{0},\mathrm{1}\right]$.

4.3 Experimental results

Our first case study involves irregularly scattered points, commonly observed in borehole datasets, with a fault dipping to the west (Fig. 5). To better understand the impact of preferred edge orientation and the imbalance between points on the hanging wall and footwall, we also examined regularly scattered data with a fault trending obliquely to the main grid axes (Fig. 6). In these case studies, the true (ground truth) orientations of the faults were 270 and 340.43°, respectively, which can be compared against the calculated statistics (Table 1).

Standard statistical methods, including the estimation of confidence intervals, demonstrated promising utility in inferring the true dip direction. For example, the minor deviations of 1–2° between error-free and with-error scenarios suggest that the model is robust in the presence of elevation errors, reinforcing the conclusions from Proposition 4. These minor deviations indicate that the method provides reliable dip direction estimates even under uncertain elevation conditions.

However, limitations remain. The confidence intervals did not contain the true dip direction, with a mismatch ranging from 2 to 5°. It is unclear whether this mismatch is due to the dataset limitations or the elevation errors' specific characteristics. Although both the directional dispersion and circular standard error increased slightly in the case studies with uncertain elevations (Table 1), the increase was not significant, and the circular standard error remained low, resulting in very narrow confidence intervals.

Figure 5 A geological horizon displaced by a vertical fault; its normal vector points to the west (azimuth 270°). (a) A case study with elevation data without errors. (b) A case study with points having elevation data with errors (mean = 0, standard deviation = 0.05). (c) A spherical projection for the orientation measurements as three-element subsets (n=122 304) of the points corresponding to panel (a). (d) A spherical projection for the orientation measurements as three-element subsets (n=161 700) of the points corresponding to panel (b). The small cloud at the center of the plot (pink polygon) corresponds to almost-flat triangles lying entirely on the same side of the fault. The spherical projection was conducted using Dips software (Rocscience, 2017). We used pole vectors, upper-hemisphere projection, and equal-angle projection.

Figure 6 Orientation measurements using a combinatorial algorithm for a regular grid of points and a vertical fault trending obliquely to the main axes of the grid. The vertical fault has a normal vector with azimuth 340.43°. (a) Input points without elevation errors. (b) Input points with elevation errors (mean = 0, standard deviation = 0.05). (c) A spherical projection for the orientation measurements as three-element subsets (n=109 111) of the points without elevation errors corresponding to panel (a). (d) A spherical projection for the orientation measurements as three-element subsets (n=157 252) of the points with elevation errors corresponding to panel (b). The small cloud at the center of the plot (pink polygon) corresponds to almost-flat triangles lying entirely on the same side of the fault. The spherical projection was conducted using Dips software (Rocscience, 2017). We used pole vectors, upper-hemisphere projection, and equal-angle projection.

Table 1 Statistical results for the case studies presented in Figs. 5 and 6. The ground truth values are 270 and 340.43° for Figs. 5 and 6, respectively.

5 Discussion

5.1 New insights from geometric and data analysis

In planar structures such as triangles, the strike is defined by the intersection of the triangle and a horizontal plane (Fossen, 2006). Therefore, a triangle with a flat-lying edge can be interpreted as having a strike parallel to this edge. The determination of the dip direction for such triangles requires solving a computational geometry problem: determining whether the third point lies to the left or right of this flat-lying edge (De Berg et al., 2008).

Our study shows that flat-lying edges parallel to the fault strike are privileged due to the greater number of points lying on one side of the edge. As a result, these directions carry more weight in statistical calculations and are more likely to represent the true dip direction. Identifying these privileged dip directions is essential for accurate fault orientation predictions, as the concentration of observations around the true dip direction indicates a reliable methodology. Previous work identified an issue of the spatial distribution of points in relation to the boundary of the study area and the fault strike (Michalak et al., 2021). However, the work here addressed this problem (Fig. 6a, b). The formal reasoning and experimental results suggests that, even if the fault strike trends obliquely to the boundary of the study area and the alignment of the regular grid of points, the algorithm still gives a preference for the edges which are parallel to the fault strike.

At this stage, it may be beneficial to complement standard statistical approaches with qualitative observations based on the distribution of dip directions. For example, in Figs. 5c, d and 6c, d, the directions opposite to the true dip direction (90 and 160.43°) are sparsely represented, which may help constrain the interpretation by highlighting directions that are unlikely to represent the true dip. This pattern aligns with the formal results from Proposition 1, which indicate that, for a fixed edge parallel to the fault strike on the upper wall (hanging wall), no dip directions will point precisely toward the upper wall (Observation 1). This observation provides additional confidence in the reliability of the predicted dip directions.

5.2 Comparison with similar approaches

In most studies related to geological model uncertainty, multiple faults are considered to estimate the uncertainty of the fault network or geometry (Cherpeau and Caumon, 2015; Lecour et al., 2001). The variability of faults is typically expressed through their parametrization, which often includes the strike and dip of the fault (Cherpeau et al., 2012; Aydin and Caers, 2017; Goodwin et al., 2022). In contrast, in this study, we set the orientation of the fault constant to provide ground truth data, enabling us to investigate the mathematical relationships between points and dip directions of triangles genetically related to the fault.

By keeping the fault orientation constant, we could isolate the effects of data uncertainties and mathematically investigate the relationship between points and dip directions. This controlled environment provides a more precise assessment of how data uncertainties, specifically elevation errors, influence dip direction calculations without the added complexity of variable fault parameters. One important limitation of the study is that the proposed method does not differentiate between normal and reverse dip-slip faults (see a broader discussion in Michalak et al., 2021).

To investigate the impact of data uncertainties on the calculated fault orientation, we added errors to the elevation values (Z coordinate), following a normal probability distribution with a mean equal to the measured elevation. The assumption of a constant Z coordinate is conceptually similar to the simple kriging approach, in which the expected value of the random variable is assumed to be the same at every point x in the domain (Wackernagel, 1995). Since our model assumes independence of errors across boreholes, a possible direction for future development would be to explore the effects of spatially correlated errors. In our case, the geographical coordinates (X and Y) are assumed to be known. Introducing uncertainties in the geographical coordinates, where locations follow a normal probability distribution (Allmendinger et al., 2011), would present a significant challenge; when coordinates are uncertain, it becomes difficult to define whether a point lies to the left or right of a line, as the vertices of that line are also uncertain.

Figure 7 Triangulated network of example bathymetric GEBCO data stored in a quasi-regular grid. (a) A repeating valley–ridge pattern with a marked rectangle analyzed in panel (b); it is located at a curved axis of a negative landform as testified by opposing dip directions near the axis. (b) A zoom-in of the rectangle marked in panel (a); here we can see that triangles with the same elevations of nodes yield a flat triangle, that triangles with a flat-lying edge but a third vertex with a different elevation than the other two points have dip direction orthogonal to the flat edge, and that the interiors of triangles are filled with the color corresponding to the value of the azimuth from panel (a). (c) A histogram of circular data related to the TIN model. The histogram shows narrow azimuthal groups with approximately 4° spacing corresponding to the edges of a regular grid. The length and color of each segment represent the number of triangles in the corresponding group.

5.3 Importance of the results in modeling spatial data

The results of this study can be applied in GIS-based directional and statistical analyses of topographic vector data, such as Triangulated Irregular Network (TIN) objects frequently used in GIS software. In particular, our findings support the explanation of singular statistical effects in azimuthal analyses of regular topographic data. To support this claim, we note the following:

1. TIN models are used in the analysis of bathymetric GEBCO (Gridded Bathymetry Data, 2024) data (Fig. 7a) (Włodarczyk-Sielicka et al., 2022; Idzikowska et al., 2024).

2. GEBCO data points are often stored in quasi-regular point grids. Although these points are regularly spaced in geographic coordinates, when projected into a Cartesian system, the individual points are no longer regularly aligned. This leads to the formation of quasi-parallel edges in the TIN model (Fig. 7b). This quasi-regular alignment of points often results in triangles with nearly parallel edges, creating preferred directions for dip calculations.

3. Limited precision of the elevation measurements can lead to the rounding of elevation values to integers. If points within a specified neighborhood are recorded as having only two distinct integer values, then this results in a constant elevation difference between two sets of points. This rounding can lead to flat edges (edges with two identical elevation values) or flat triangles (triangles with three identical elevation values), as marked in Fig. 7b.

For the specific combination of models, spatial distribution of points, and limited precision described above, the conditions from Proposition 1 may apply, explaining the concentration of dip direction values in azimuthal histograms (Fig. 7c). Only eight narrow groups of azimuth values are present in the entire dataset due to the eight preferred triangle edge directions in the TIN. Two groups (130 and 313°) are significantly underrepresented, likely due to the lack of flat edges striking NE–SW. This feature is purely based on the data representation. Users must take this into account during analysis to avoid skewed interpretations of fault orientations.

6 Conclusion

This study aimed to bridge computational geometry and structural geology to explore the behavior of triangles related to a horizon displaced by a vertical fault. We conducted analyses under both idealized (two elevations) and uncertain (added elevation error) conditions. A key challenge was adapting concepts from computational geometry, such as determining whether a point lies to the left or right of a line, to a geological context. The main conclusion is particularly useful for researchers working with triangulated surfaces – either explicitly or implicitly – for directional surface analysis (e.g., azimuth, dip direction maps), especially when those surfaces are affected by sampling uncertainty.

Key findings from the study include the following:

• This study shows that, for directional analyses of triangulated surfaces, determining the dip direction of a triangle with a flat-lying edge requires us to solve a computational geometry problem, determining whether the third point lies to the left or right of the edge.

• The problem of assessing fault orientation can be approached as an optimization task, where the fault orientation is estimated by identifying the edge with the maximum number of points on one side, leading to the maximum number of triangles with the same dip direction.

• Our formal analysis shows that introducing measurement errors does not affect the expected dip direction of samples, which remains identical to the error-free case. Moreover, the statistical results and orientation distributions remain robust across different fault and study area orientations, suggesting practical applications.

• The importance of these findings is particularly relevant for directional analyses of imprecise topographic data, such as azimuth maps of bathymetric point datasets distributed in regular grids. The concentration of azimuths around the N–S, W–E, NE–SE, and NE–SW directions corresponds directly to the edges of the triangulation in the regular grid.

Appendix A

To generate all possible triangles from a given set of boreholes, we used an algorithm to generate an all k-element (k=3) subsets from an n-element set X (k<n) in a lexicographic order (Lipski, 2004). To explain how the algorithm works, we first note that every k-element subset can be uniquely represented by an increasing sequence of length k of elements from X. For example, the subset $\mathit{\{}\mathrm{3},\mathrm{5},\mathrm{1}\mathit{\}}$ can be represented as a sequence (1, 3, 5). The first step in the algorithm involves writing the first k digits from X. For example, if k=4, the first sequence would be (1, 2, 3, 4). Then, the sequence succeeding $\left(a_{\mathrm{1}},\mathrm{\dots },a_k\right)$ is

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\left(b_{\mathrm{1}},\mathrm{\dots },b_k\right)& =\left(a_{\mathrm{1}},\mathrm{\dots },a_{p-\mathrm{1}},a_p+\mathrm{1},a_p+\mathrm{2},\right.\\ & \left.\mathrm{\dots },a_p+k-p+\mathrm{1}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A2)}& p=max\mathit{\{}i:a_i<n-k+i\mathit{\}}.\end{array}$$

Likewise, the sequence which succeeds $\left(b_{\mathrm{1}},\mathrm{\dots },b_k\right)$ is

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\left(c_{\mathrm{1}},\mathrm{\dots },c_k\right)& =\left(b_{\mathrm{1}},\mathrm{\dots },b_{p^{\prime }-\mathrm{1}},b_{p^{\prime }}+\mathrm{1},b_{p^{\prime }}+\mathrm{2},\mathrm{\dots },b_{p^{\prime }}\right.\\ & \left.+k-p^{\prime }+\mathrm{1}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& p^{\prime }=p-\mathrm{1}\text{ when }{b}_k=n\text{ and}\\ & p^{\prime }=k\text{ when }{b}_k<n.\end{array}\end{array}$$

During the procedure, we assume that the sequences $\left(a_{\mathrm{1}},\mathrm{\dots },a_k\right)$ and $\left(b_{\mathrm{1}},\mathrm{\dots },b_k\right)$ are different from $\left(n-k+\mathrm{1},\mathrm{\dots },n\right)$, the last sequence in our order. Here, p and p^′ can be conceptualized as indices where updates of digits starting from the largest (k) index terminate. If p or p^′ are determined, then a_p or $b_{p^{\prime }}$ can easily be found, allowing their use in the update procedure (Eqs. A1 and A3). We note that the number of k-element subsets from an n-element set can be determined using the binomial coefficient. For example, if n=6 and k=4, then

$$\begin{array}{}\text{(A5)}& \left(\begin{array}{l}n\\ k\end{array}\right)={\displaystyle \frac{n\mathrm{!}}{k\mathrm{!}\left(n-k\right)\mathrm{!}}}={\displaystyle \frac{\mathrm{4}\mathrm{!}\cdot \mathrm{5}\cdot \mathrm{6}}{\mathrm{4}\mathrm{!}\cdot \mathrm{2}\mathrm{!}}}={\displaystyle \frac{\mathrm{30}}{\mathrm{2}}}=\mathrm{15}.\end{array}$$

Appendix B: A formal proof of Proposition 1

Proposition 1.

Proof. We note that, intuitively, the result from parts A and B is in line with the standard concept of a strike of a geological planar feature defined as the intersection of this planar feature and a horizontal plane (Fossen, 2006; see also Fig. 4 as an illustration of the proposition).

To formally prove the proposition, particularly part C, we will firstly refer to the following orientation test (De Berg et al., 2008).

Fact 1 (orientation test) (De Berg et al., 2008).

Given three vertices, $s=\left(s_{\mathrm{1}},s_{\mathrm{2}}\right),t=\left(t_{\mathrm{1}},t_{\mathrm{2}}\right)$, and $u=\left(u_{\mathrm{1}},u_{\mathrm{2}}\right)$, the sign of the determinant

$$\begin{array}{}\text{(B1)}& \left|\mathbf{D}\right|=\text{det}\left[\begin{array}{rrr}\mathrm{1}& t_{\mathrm{1}}& t_{\mathrm{2}}\\ \mathrm{1}& u_{\mathrm{1}}& u_{\mathrm{2}}\\ \mathrm{1}& s_{\mathrm{1}}& s_{\mathrm{2}}\end{array}\right]\end{array}$$

determines whether s lies left or right of the oriented vector tu.

Proof of Proposition 1.

Part A of Proposition 1.

Let $p_{\mathrm{1}}=\left(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}\right)$ and $p_{\mathrm{2}}=\left(x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}}\right)$ be points forming an edge. For simplicity, let us assume that p₁ and p₂ are located on the surface of the footwall. The third point $p_{\mathrm{3}}=\left(x_{\mathrm{3}},y_{\mathrm{3}},z_{\mathrm{3}}\right)$ can be anywhere on the surface of the hanging wall. Therefore, we will consider a set of many triangles. Because the walls are horizontal and the difference between walls is constant, we can write $p_{\mathrm{1}}=\left(x_{\mathrm{1}},y_{\mathrm{1}},a\right)$, $p_{\mathrm{2}}=\left(x_{\mathrm{2}},y_{\mathrm{2}},a\right)$, and $p_{\mathrm{3}}=\left(x_{\mathrm{3}},y_{\mathrm{3}},b\right)$, where a≠b. Let $k:=a-b$ be the positive constant, the elevation difference between walls.

The two vectors spanning the triangle's plane are as follows:

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}},\mathrm{0}\right],\text{(B3)}& {\displaystyle }{\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{3}}-x_{\mathrm{1}},y_{\mathrm{3}}-y_{\mathrm{1}},-k\right].\end{array}$$

Using the cross-product, the coordinates of the normal vector of the triangle defined by the above vectors can be calculated via the mnemonic rule:

$$\begin{array}{}\text{(B4)}& \mathrm{\text{det}}\left[\begin{array}{lll}{x}_{\mathrm{2}}-x_{\mathrm{1}}& x_{\mathrm{3}}-x_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{1}\right]\\ y_{\mathrm{2}}-y_{\mathrm{1}}& y_{\mathrm{3}}-y_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{2}\right]\\ \mathrm{0}& -k& n_{\mathrm{1}}\left[\mathrm{3}\right]\end{array}\right]=\begin{array}{l}\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\\ \cdot n_{\mathrm{1}}\left[\mathrm{3}\right]+\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ \cdot -k\cdot n_{\mathrm{1}}\left[\mathrm{1}\right]+\mathrm{0}\\ \cdot (x_{\mathrm{3}}-x_{\mathrm{1}})\cdot n_{\mathrm{1}}\left[\mathrm{2}\right]\\ -n_{\mathrm{1}}\left[\mathrm{1}\right]\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\\ \cdot \mathrm{0}-n_{\mathrm{1}}\left[\mathrm{2}\right]\cdot -k\\ \cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)-n_{\mathrm{1}}\left[\mathrm{3}\right]\\ \cdot (x_{\mathrm{3}}-x_{\mathrm{1}})\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}\end{array}$$

In summary,

$$\begin{array}{}\text{(B5)}& {\displaystyle }\begin{array}{rl}{n}_{\mathrm{1}}\left[\mathrm{1}\right]& =\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\cdot -k-\left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\cdot \mathrm{0}\\ & =\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\cdot -k\\ & =-k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\end{array}\text{(B6)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{2}\right]=\mathrm{0}\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)-\left(-k\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)=k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\text{(B7)}& {\displaystyle }\begin{array}{rl}{n}_{\mathrm{1}}\left[\mathrm{3}\right]& =\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)-\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ & =x_{\mathrm{2}}{y}_{\mathrm{3}}-x_{\mathrm{2}}{y}_{\mathrm{1}}-x_{\mathrm{1}}{y}_{\mathrm{3}}+x_{\mathrm{1}}{y}_{\mathrm{1}}-x_{\mathrm{3}}{y}_{\mathrm{2}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}\\ & =x_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}-x_{\mathrm{2}}{y}_{\mathrm{1}}-x_{\mathrm{1}}{y}_{\mathrm{3}}-x_{\mathrm{3}}{y}_{\mathrm{2}}.\end{array}\end{array}$$

As of now, we have two notes:

• Every non-vertical triangle has two normal vectors: one is directed downwards, and the other is directed upwards. We are only interested in normal vectors directed upwards to avoid duplicate representations and ensure consistent representation of observations.

• Note that only the edge on the footwall is fixed and not the third point, whose coordinates affect the sign of the third coordinate of the normal vector (Eq. B7). Therefore, it can be concluded that we do not investigate a particular normal vector but a set of many normal vectors. Moreover, given only the position of the fixed edge, the normal vector's third coordinate is unknown, and three cases must be considered: when it is positive, negative, or zero. We will now consider these three cases.

• I.

  If n₁[3]>0, then the coordinates of the normal vector look as above, which means that the vector is directed upwards.

• II.

  If n₁[3]<0, then it means that the normal vector is directed downwards, and the coordinates must be multiplied by minus 1 to adhere to the above rule that there is only one representation of normal vectors.

• III.

  If n₁[3]=0, then the normal vector is directed horizontally (orthogonal to the vectors ${\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}}\right]$ and ${\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{3}}-x_{\mathrm{1}},y_{\mathrm{3}}-y_{\mathrm{1}}\right]$), and the corresponding triangle is vertical, contrary to the initial assumption of considering only non-vertical triangles. Therefore, we no longer consider this scenario.

Now, according to (II), we multiply coordinates from Eqs. (B5)–(B7) by minus 1.

$$\begin{array}{}\text{(B8)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{1}\right]=k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\text{(B9)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{2}\right]=-k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\text{(B10)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{3}\right]=-\left(x_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}-x_{\mathrm{2}}{y}_{\mathrm{1}}-x_{\mathrm{1}}{y}_{\mathrm{3}}-x_{\mathrm{3}}{y}_{\mathrm{2}}\right)\end{array}$$

We observe that the coordinates of normal vectors (Eqs. B5–B7 and B8–B10) are not the same, which means that we obtained two distinct normal vectors directed upwards, n₁ and n₂, and which have a dip direction difference of 180°. The vectors are directed upwards because the third coordinate is positive. So, part A of Proposition 1 is proven.

Part B of Proposition 1.

As we already know from (A), there are only two possible dip directions of the infinite set of triangles. These dip directions have a directional difference of 180°. Therefore, to prove that they are orthogonal to the edge ${\mathbf{p}}_{\mathrm{1}}{\mathbf{p}}_{\mathrm{2}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}}\right]$, it is enough to prove this for one vector: n₁ or n₂ projected on the horizontal plane (the projection of the normal vector onto the horizontal plane does not change its direction). Using the dot product (⋅), we can show it for ${\stackrel{\mathrm{˘}}{\mathit{n}}}_{\mathrm{1}}$, where ${\stackrel{\mathrm{˘}}{\mathit{n}}}_{\mathrm{1}}$ denotes the projection (a 2D variant) of the first normal vector onto the horizontal plane:

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{\mathbf{p}}_{\mathrm{1}}{\mathbf{p}}_{\mathrm{2}}\cdot {\stackrel{\mathrm{˘}}{\mathit{n}}}_{\mathrm{1}}& =\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot -k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)+(y_{\mathrm{2}}-y_{\mathrm{1}})\cdot k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\\ & =-k\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)+k\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)=\mathrm{0}.\end{array}\end{array}$$

Part C of Proposition 1.

We have a fixed edge on the footwall, and we can consider two cases: when the point lies to the left or to the right of the line containing the edge.

Using Fact 1 (orientation test), we can determine whether $p_{\mathrm{3}}=\left(x_{\mathrm{3}},y_{\mathrm{3}}\right)$ lies to the left, to the right, or on the edge p₁p₂:

$$\begin{array}{}\text{(B12)}& \left|\mathbf{D}\right|=\text{det}\left[\begin{array}{rll}\mathrm{1}& x_{\mathrm{1}}& y_{\mathrm{1}}\\ \mathrm{1}& x_{\mathrm{2}}& y_{\mathrm{2}}\\ \mathrm{1}& x_{\mathrm{3}}& y_{\mathrm{3}}\end{array}\right]=\begin{array}{l}{x}_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}\\ -y_{\mathrm{1}}{x}_{\mathrm{2}}-y_{\mathrm{2}}{x}_{\mathrm{3}}-y_{\mathrm{3}}{x}_{\mathrm{1}}.\end{array}\end{array}$$

We note that the value can be positive, negative, or zero depending on whether the point p₃ lies left, right, or on the edge p₁p₂. Recall that, in part A, this value was also the third coordinate of the first non-horizontal normal vector n₁[3], and the additive inverse $-\left(x_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}-x_{\mathrm{2}}{y}_{\mathrm{1}}-x_{\mathrm{1}}{y}_{\mathrm{3}}-x_{\mathrm{3}}{y}_{\mathrm{2}}\right)$ was the third coordinate of the second normal vector. Therefore, the signs of the expressions $\text{sgn}\left(x_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}-y_{\mathrm{1}}{x}_{\mathrm{2}}-y_{\mathrm{2}}{x}_{\mathrm{3}}-y_{\mathrm{3}}{x}_{\mathrm{1}}\right)\ne \text{sgn}\left(-\left(x_{\mathrm{2}}{y}_{\mathrm{3}}+x_{\mathrm{3}}{y}_{\mathrm{1}}+x_{\mathrm{1}}{y}_{\mathrm{2}}-y_{\mathrm{1}}{x}_{\mathrm{2}}-y_{\mathrm{2}}{x}_{\mathrm{3}}-y_{\mathrm{3}}{x}_{\mathrm{1}}\right)\right)$ simultaneously determine the position of the point p₃ relative to the edge p₁p₂ and the choice of one of two possible normal vectors. In the case of $\left|D\right|=\mathrm{0}$, the points p₁,p₂ (footwall) and p₃ (hanging wall) are collinear in 2D space but not in 3D space, and the corresponding triangle is vertical, which is beyond the scope of our study.

Appendix C

Proposition 3. Analysis with elevation uncertainties restricted to the hanging wall regarding the free vertices.

Similar calculations can be conducted for data with fixed geographical position but uncertain elevations. To achieve this, the uncertain elevations can be represented as sums of the measured constant elevations and a random variable ε with the normal distribution N(0,σ²).

Therefore, the uncertain elevations of n points are independent random variables ${\mathit{\epsilon }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\epsilon }}_n$, and their expected values are equal to zero; i.e., $E\left[{\mathit{\epsilon }}_{\mathrm{1}}\right]=\mathrm{\dots }=E\left[{\mathit{\epsilon }}_n\right]=\mathrm{0}$. Here we consider $x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}},x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}}$ to be fixed constants.

From a practical viewpoint, we firstly create the set of points from the uniform distribution displaced by a vertical fault and then add error to elevation data.

From then on, we have points $p_{\mathrm{1}}=\left(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{1}}\right)$ and $p_{\mathrm{2}}=\left(x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{2}}\right)$ forming an edge on the footwall. The third point $p_{\mathrm{3}}=\left(x_{\mathrm{3}},y_{\mathrm{3}},z_{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{3}}\right)$ can lie anywhere on the hanging wall. For now, we consider the point p₃ to be fixed, but, ultimately, it will traverse the points on the hanging wall to compute their average.

The random vectors spanning the plane of a random triangle are calculated as follows:

$$\begin{array}{}\text{(C1)}& {\displaystyle }\begin{array}{rl}& {\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}}-{\mathit{\epsilon }}_{\mathrm{1}}\right],\text{ and let }\\ & \mathit{\alpha }:={\mathit{\epsilon }}_{\mathrm{2}}-{\mathit{\epsilon }}_{\mathrm{1}}.\text{ Then }\\ & {\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}},\mathit{\alpha }\right],\end{array}\text{(C2)}& {\displaystyle }\begin{array}{rl}& {\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{3}}-x_{\mathrm{1}},y_{\mathrm{3}}-y_{\mathrm{1}},z_{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{3}}-z_{\mathrm{1}}-{\mathit{\epsilon }}_{\mathrm{1}}\right],\text{ and let }\\ & k:=-\left(z_{\mathrm{2}}-z_{\mathrm{1}}\right)\text{ with }\mathit{\beta }:={\mathit{\epsilon }}_{\mathrm{3}}-{\mathit{\epsilon }}_{\mathrm{1}}.\\ & \text{Then }{\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{3}}-x_{\mathrm{1}},y_{\mathrm{3}}-y_{\mathrm{1}},-k+\mathit{\beta }\right].\end{array}\end{array}$$

The normal vector can be calculated using the mnemonic rule for cross-product:

$$\begin{array}{}\text{(C3)}& \mathrm{\text{det}}\left[\begin{array}{lll}{x}_{\mathrm{2}}-x_{\mathrm{1}}& x_{\mathrm{3}}-x_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{1}\right]\\ y_{\mathrm{2}}-y_{\mathrm{1}}& y_{\mathrm{3}}-y_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{2}\right]\\ \mathit{\alpha }& -k+\mathit{\beta }& n_{\mathrm{1}}\left[\mathrm{3}\right]\end{array}\right]=\begin{array}{l}\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\\ \cdot n_{\mathrm{1}}\left[\mathrm{3}\right]+\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ \cdot \left(-k+\mathit{\beta }\right)\cdot n_{\mathrm{1}}\left[\mathrm{1}\right]\\ +\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot n_{\mathrm{1}}\left[\mathrm{2}\right]\\ -n_{\mathrm{1}}\left[\mathrm{1}\right]\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\cdot \mathit{\alpha }\\ -\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(-k+\mathit{\beta }\right)\\ \cdot n_{\mathrm{1}}\left[\mathrm{2}\right]-\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ \cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot n_{\mathrm{1}}\left[\mathrm{3}\right].\end{array}\end{array}$$

In summary,

$$\begin{array}{}\text{(C4)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{1}\right]=\left(-k+\mathit{\beta }\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\text{(C5)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{2}\right]=\left(k-\mathit{\beta }\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\text{(C6)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{3}\right]=\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)-\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}$$

As in the deterministic case, we note that the third coordinate can be positive or negative. If it is positive, then the coordinates stay as they are. Otherwise, all coordinates must be multiplied by minus 1. Then, the coordinates are as follows:

$$\begin{array}{}\text{(C7)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{1}\right]=\left(k-\mathit{\beta }\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\text{(C8)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{2}\right]=\left(-k+\mathit{\beta }\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\text{(C9)}& {\displaystyle }{n}_{\mathrm{2}}\left[\mathrm{3}\right]=-\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)+\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}$$

For both cases, we can now calculate the expectations of the first two coordinates, given that the above expressions contain random variables.

For the first case (positive Z value),

$$\begin{array}{}\text{(C10)}& \begin{array}{l}E\left[n_{\mathrm{1}}\left[\mathrm{1}\right]\right]=\\ E\left[\left(-k+\mathit{\beta }\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\right]=\\ (-k+E\left[\mathit{\beta }\left]\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-E[\mathit{\alpha }\right]\cdot (y_{\mathrm{3}}-y_{\mathrm{1}})=\\ \left(-k+\left(E\left[{\mathit{\epsilon }}_{\mathrm{3}}]-E[{\mathit{\epsilon }}_{\mathrm{1}}\right]\right)\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ -\left(E\left[{\mathit{\epsilon }}_{\mathrm{2}}]-E[{\mathit{\epsilon }}_{\mathrm{1}}\right]\right)\cdot (y_{\mathrm{3}}-y_{\mathrm{1}})=\\ -k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}\end{array}$$

Since E[n₁[1]] does not depend on y₃, averaging over all hanging wall points does not change the expected value, which will remain equal to $-k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)$:

$$\begin{array}{}\text{(C11)}& \begin{array}{l}E\left[n_{\mathrm{1}}\left[\mathrm{2}\right]\right]=\\ E\left[\left(k-\mathit{\beta }\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\right]=\\ \left(k-\left(E\left[{\mathit{\epsilon }}_{\mathrm{3}}]-E[{\mathit{\epsilon }}_{\mathrm{1}}\right]\right)\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\\ +\left(E\left[{\mathit{\epsilon }}_{\mathrm{2}}]-E[{\mathit{\epsilon }}_{\mathrm{1}}\right]\right)\cdot (x_{\mathrm{3}}-x_{\mathrm{1}})=\\ k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right).\end{array}\end{array}$$

Since E[n₁[2]] does not depend on x₃, averaging over all hanging wall points does not change the expected value, which will remain equal to $k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)$.

For the second case (negative Z value),

$$\begin{array}{}\text{(C12)}& \begin{array}{l}E\left[n_{\mathrm{2}}\left[\mathrm{1}\right]\right]=\\ E\left[\left(k-\mathit{\beta }\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\right]=\\ k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}\end{array}$$

As previously, since E[n₂[1]] does not depend on y₃, averaging over all hanging wall points does not change the expected value, which will remain equal to $k\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)$.

$$\begin{array}{}\text{(C13)}& \begin{array}{l}E\left[n_{\mathrm{2}}\left[\mathrm{2}\right]\right]=\\ E\left[\left(-k+\mathit{\beta }\right)\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\right]=\\ -k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\end{array}\end{array}$$

Again, as E[n₂[2]] does not depend on x₃, averaging over all hanging wall points does not change the expected value, which will remain equal to $-k\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)$.

Appendix D

Proposition 4. Analysis with elevation uncertainties with free points on the surfaces of the hanging wall or footwall.

Here, we continue the analysis from the Sect. 4.2 by considering only triangles with third points lying on the surface of the footwall (from v_k+1 to v_n). As we have n−k triangles, we need $n-k+$2 points on the footwall. Denote them by $p_i=\left(x_i,y_i,z_i+{\mathit{\epsilon }}_i\right)$, where $i=\mathrm{1},\mathrm{2},\mathrm{\dots },n-k+\mathrm{2}$.

We fix the edge e on the surface of the footwall defined by two points: $p_{\mathrm{1}}=\left(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{1}}\right)$ and $p_{\mathrm{2}}=\left(x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{2}}\right)$. Because the edge e is fixed, the associated coordinates $x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}},x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}}$ are considered fixed constants. The third point p₃ is sampled from the surface of the footwall. As previously, we fix the point p₃ for a moment, but, ultimately, it will traverse the points on the hanging wall to compute their average.

The third point $p_{\mathrm{3}}=\left(x_{\mathrm{3}},y_{\mathrm{3}},z_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{3}}\right)$ is sampled from the surface of the footwall. We calculate the spanning vectors

$${\mathit{v}}_{\mathrm{1}}=\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}},\mathit{\alpha }\right]$$

and

$${\mathit{v}}_{\mathrm{2}}=\left[x_{\mathrm{3}}-x_{\mathrm{1}},y_{\mathrm{3}}-y_{\mathrm{1}},\mathit{\beta }\right],$$

with

$$\mathit{\alpha }:={\mathit{\epsilon }}_{\mathrm{2}}-{\mathit{\epsilon }}_{\mathrm{1}}$$

and

$$\mathit{\beta }:={\mathit{\epsilon }}_{\mathrm{3}}-{\mathit{\epsilon }}_{\mathrm{1}}.$$

The normal vector can be calculated using the mnemonic rule for cross-product:

$$\begin{array}{}\text{(D1)}& \mathrm{\text{det}}\left[\begin{array}{lll}{x}_{\mathrm{2}}-x_{\mathrm{1}}& x_{\mathrm{3}}-x_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{1}\right]\\ y_{\mathrm{2}}-y_{\mathrm{1}}& y_{\mathrm{3}}-y_{\mathrm{1}}& n_{\mathrm{1}}\left[\mathrm{2}\right]\\ \mathit{\alpha }& \mathit{\beta }& n_{\mathrm{1}}\left[\mathrm{3}\right]\end{array}\right]=\begin{array}{l}\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\\ \cdot n_{\mathrm{1}}\left[\mathrm{3}\right]+\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ \cdot \left(\mathit{\beta }\right)\cdot n_{\mathrm{1}}\left[\mathrm{1}\right]+\mathit{\alpha }\\ \cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot n_{\mathrm{1}}\left[\mathrm{2}\right]\\ -n_{\mathrm{1}}\left[\mathrm{1}\right]\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\\ \cdot \mathit{\alpha }-\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(\mathit{\beta }\right)\\ \cdot n_{\mathrm{1}}\left[\mathrm{2}\right]-\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)\\ \cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot n_{\mathrm{1}}\left[\mathrm{3}\right].\end{array}\end{array}$$

In summary,

$$\begin{array}{}\text{(D2)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{1}\right]=\mathit{\beta }\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\text{(D3)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{2}\right]=-\mathit{\beta }\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\text{(D4)}& {\displaystyle }{n}_{\mathrm{1}}\left[\mathrm{3}\right]=\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)-\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}$$

As previously, we start the calculations of expected values of E[n₁[1]], E[n₁[2]], E[n₁[3]] and then the averages over all footwall points.

$$\begin{array}{}\text{(D5)}& \begin{array}{l}E\left[n_{\mathrm{1}}\left[\mathrm{1}\right]\right]=\\ E\left[\mathit{\beta }\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-\mathit{\alpha }\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)\right]=\\ E\left[\mathit{\beta }\right]\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)-E\left[\mathit{\alpha }\right]\left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)=\mathrm{0}\end{array}\end{array}$$

Since E[n₁[1]]=0, the average over all footwall points does not change and will remain equal to 0.

$$\begin{array}{}\text{(D6)}& \begin{array}{l}E\left[n_{\mathrm{1}}\left[\mathrm{2}\right]\right]=\\ E\left[-\mathit{\beta }\cdot \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+\mathit{\alpha }\cdot \left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\right]=\\ -E\left[\mathit{\beta }\right]\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+E\left[\mathit{\alpha }\right]\left(E\left[x_{\mathrm{3}}-x_{\mathrm{1}}\right]\right)=\mathrm{0}\end{array}\end{array}$$

Since E[n₁[2]]=0, the average over all footwall points does not change and will remain equal to 0.

$$\begin{array}{}\text{(D7)}& \begin{array}{l}E\left[n_{\mathrm{1}}\left[\mathrm{3}\right]\right]=\\ \left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)-\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)\cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right).\end{array}\end{array}$$

Now, taking the average over all footwall points, we get

$$\begin{array}{}\text{(D8)}& \begin{array}{rl}\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)& \cdot {\displaystyle \frac{\left(y_{\mathrm{3}}-y_{\mathrm{1}}\right)+\mathrm{\dots }+\left(y_{n-k+\mathrm{2}}-y_{\mathrm{1}}\right)}{n-k}}\\ & -{\displaystyle \frac{(x_{\mathrm{3}}-x_{\mathrm{1}})+\mathrm{\dots }+(x_{n-k+\mathrm{2}}-x_{\mathrm{1}})}{n-k}}\\ & \cdot (y_{\mathrm{2}}-y_{\mathrm{1}}).\end{array}\end{array}$$

Observe that E[n₁[3]]=0 if and only if

$$\begin{array}{}\text{(D9)}& \begin{array}{rl}\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)& \cdot \left[\right(y_{\mathrm{3}}-y_{\mathrm{1}})+\mathrm{\dots }+(y_{n-k+\mathrm{2}}-y_{\mathrm{1}}\left)\right]\\ & =\left[\right(x_{\mathrm{3}}-x_{\mathrm{1}})+\mathrm{\dots }+(x_{n-k+\mathrm{2}}-x_{\mathrm{1}}\left)\right]\\ & \cdot \left(y_{\mathrm{2}}-y_{\mathrm{1}}\right),\end{array}\end{array}$$

i.e., if and only if either $\mathbf{w}=\left[\right(x_{\mathrm{3}}-x_{\mathrm{1}})+\mathrm{\dots }+(x_{n-k+\mathrm{2}}-x_{\mathrm{1}})$, $(y_{\mathrm{3}}-y_{\mathrm{1}})+\mathrm{\dots }+(y_{n-k+\mathrm{2}}-y_{\mathrm{1}})]$ is the zero vector or it is non-zero and parallel to the vector $\left[x_{\mathrm{2}}-x_{\mathrm{1}},y_{\mathrm{2}}-y_{\mathrm{1}}\right]$, which happens very seldom.

One of the trivial cases for which E[n₁[3]]≠0 is when

$$\begin{array}{}\text{(D10)}& \begin{array}{rl}& {\displaystyle \frac{\left(x_{\mathrm{3}}-x_{\mathrm{1}}\right)+\mathrm{\dots }+(x_{n-k+\mathrm{2}}-x_{\mathrm{1}})}{x_{\mathrm{2}}-x_{\mathrm{1}}}}\\ & \ne {\displaystyle \frac{(y_{\mathrm{3}}-y_{\mathrm{1}})+\mathrm{\dots }+(y_{n-k+\mathrm{2}}-y_{\mathrm{1}})}{y_{\mathrm{2}}-y_{\mathrm{1}}}},\end{array}\end{array}$$

assuming x₂≠x₁ and y₂≠y₁.

Another trivial case for which E[n₁[3]]≠0 is when x₂=x₁ and $(x_{\mathrm{3}}-x_{\mathrm{1}})+\mathrm{\dots }+(x_{n-k+\mathrm{2}}-x_{\mathrm{1}})\ne$ or when y₂=y₁ and $(y_{\mathrm{3}}-y_{\mathrm{1}})+\mathrm{\dots }+(y_{n-k+\mathrm{2}}-y_{\mathrm{1}})\ne \mathrm{0}$.