Mapping faults in 3D seismic data – why the method matters

3D seismic reflection imagery is the most widely used resource for interpreting the geometric structure of faults in the subsurface. Yet these endeavours carry uncertainties, the significance of which are rarely discussed. We explore how the application of different workflows yield different interpretations of a single high-quality 3D seismic image-set. We describe and apply five mapping workflows, based on 2D derivations of imagery to map an array of small-scale faults. Some workflows use vertical profiles, a conventional approach, others use plan views. We also vary the amount of under-sampling. Stacking the fault maps derived from the five workflows into a heat map shows broadly similar trends and distributions of faults irrespective of the workflow deployed. However, juxtaposition mapping reveals differences in fault length and network pattern (linkage and segmentation) arising from the different workflows. We quantify the total fault areas and distribution of fault lengths for each workflow – revealing significant differences in these statistics. Our results show that mapping strategies impact the interpretation of fault geometry, their network patterns and derived statistics. This understanding is critical for assessing and risking fault interpretations – deploying multiple workflows can reveal inherent uncertainty in structural interpretation of 3D seismic imagery.


Introduction
Geological interpretations of images derived from 3D seismic surveys provide, by far, the most extensive understanding of the structure of individual faults and the geometry of fault arrays in the subsurface.But how dependent is this understanding on the methods deployed by interpreters to map these structures?The aim of this paper is to examine how different workflows yield different 3D fault maps in a single 3D seismic volume.We concentrate on relatively small faults rather than those that bound sedimentary basins, not least because small structures have special significance for assessments of the integrity of subsurface geo-storage sites.Our work builds on some recent studies that show how under-sampling, a routine approach used to reduce the time required to work with large seismic volumes, can yield unreliable fault maps (Cunningham et al., 2021;Michie et al., 2021;Ze and Alves, 2016;2019).We confirm these results but expand the study by analysing multiple workflow approaches, alongside sampling strategy.Critical to all these approaches is the notion that the true 3D geometry of faults and fault systems in the real subsurface is always unknown, given the resolution and imaging challenges of seismic data.Therefore, beyond consideration of under-sampling, it is misleading to suggest that any particular mapping approach yields better (i.e., more accurate) 3D fault maps.

Context
In the great majority of published interpretations of fault patterns derived from seismic imagery, the workflows are undocumented, and fault maps are presented as single interpretations.In part, this deterministic approach arises from long-standing commercial drivers that underlie subsurface studies but also historical traditions of scientific advocacy (Bond, 2015;Bond et al., 2008).The exploration and development of subsurface resources such as hydrocarbons will always require the acquisition of subsurface information after the initial geological interpretations, be that from further surveys, well penetrations or production data.In practice, fault maps can be refined during this process.However, for many applications associated with the demands of the energy transition, these types of interventions and an iterative approach are not appropriate.Rather, the integrity of particular geological targets as deep storage or disposal sites (e.g., for CO 2 , H 2 or radioactive waste) require risk assessments for a broad diversity of interpretations.These demands should embrace interpretation uncertainty and therefore require the generation of arrays of different interpretations (e.g., Cherpeau et al., 2010;Cherpeau and Caumon, 2015;Schweizer et al., 2017;Godefroy et al., 2019;Han et al., 2022;Goodwin et al., 2022;Dashti et al., 2023).Our underlying motivation here is to promote not only the careful documentation of mapping workflows but also the application of different methodsand certainly not to recommend a specific workflow as being more accurate (better) than another.As the actual geological structure of the real subsurface is unknown, such designations instil overconfidence and commensurate false evaluations of site failure risks.
Uncertainty in fault maps created using seismic imagery arises from two broad areas.The first is the derivation of the seismic image from the original seismic recordings.The pathway to generate seismic reflection image datasets involves several steps with multiple alternative algorithms that give this type of imagery a non-unique nature and make their interpretation inherently uncertain (e.g., Alcalde et al., 2017;Goloshubin et al., 2006).Additionally, as with every image dataset, there are resolution constraints, spatial limitations and variable noise content that hamper interpretation.The second area of uncertainty comes from the interpretation process in which Interpreters' knowledge, prior experience and knowledge of conceptual models can play significant roles.The interpretation of faults in seismic reflection images is a challenging, decision-making process, and these decisions play a crucial role that impacts upon the final outcome (e.g., Bond, 2015;Bond et al., 2007;2008;Schaaf and Bond, 2019).As part of the interpretation process, the interpreter must make choices about the workflow they will use for interpretation, which is in part constrained and influenced by the choice of the interpretation software package.It is the aspect of workflow choice that we focus on here, considering options of view and seismic image realisation.
The interpretation of faults using seismic imagery involves a series of steps, actions, and inputs that may be selected and modified by the interpreter.Software packages offer multiple display options within which to perform fault interpretation, including a variety of tools, views, and display settings.Multiple displays give interpreters different perspectives on possible fault structures.The visualisation of 3D seismic surveys is possible in two main directions -namely via 2D slices, side-toside through vertical profiles and top-down/bottom-up using plan views (Fig. 1).It is possible to populate these views with different seismic attributes (i.e., Seismic Amplitude, Variance, Dominant frequency, etc.), and scales are adjustable (i.e., vertical exaggeration).Although seismic amplitude displays are routinely used to trace out stratal reflectors and faults, in profiles, other seismic attributes such as variance or dip are available and applied, especially in plan view.It has been shown that seismic attributes are excellent tools for highlighting discontinuities within seismic volumes (e.g., Kim et al., 2021 and references therein); however, their use requires a good understanding of attribute properties (Marfurt and Alves, 2015;Townsend et al., 1998).According to the display settings chosen by the interpreter, fault recognition and tracking criteria can vary, reinforcing the high uncertainty that the interpretation of faults carries.Interpreters also determine the number and spacing between views.The interpreter will almost always underuse the data and view options available due to time constraints often associated with industry deadlines and limited economic resources.Because interpreters can run distinct methods to interpret faults in the subsurface, we consider documentation of workflows to be especially important, and descriptions should include the interpreter's choices and criteria.The multiple options (ways) available to perform fault seismic interpretation and tight time constraints for any particular study, plus the partial lack of analysis about the impact of following determined methods, motivated this work.In the same way that the generation of single interpretations can cloud our understanding of structures in the subsurface, it is misleading to think there is a perfect/ideal method to perform fault seismic interpretation.Rather, we recommend that interpreters should provide workflows alongside their interpretations and be aware of how their choices of method can influence the resultant fault maps.

Our study
In this work, we analysed five fault mappings generated by applying different methods to the same dataset.Readers can find detailed descriptions of the methods employed in this study in section 4. Our fault interpretation methods include the employment of different types of views displaying seismic amplitude or variance attributes, distinct fault placement criteria and varied sampling density.Two interpretations use serial profiles with different spacing, one is entirely based on the use of plan views, and two combine plan views with serial profiles.We integrated the resulting five fault maps through construction in plan view of a juxtaposition map and a heat map, displaying the areas where the most fault interpretations coincide.The composite heat map is a single-composed solution which shows the locations of overlapping faults and fault patterns generated from the different workflows, which could be considered as a 'most likely' fault indicator, but we advocate consideration of the full set of mapped interpretations to give an appreciation of the uncertainty in fault mapping in the volume and to inform risk analysis.To quantitively examine the differences in fault maps derived from the workflows, we compared two fault dimensional parameters -length and area.The implications of the differences in fault patterns and parameters for each of the workflows we used, further analyses, suggestions and applications are presented in the discussion.It is not our intention to recommend any one of these workflows over the others as yielding better three-dimensional interpretations: as with any subsurface seismic mapping task, such insights are unrealistic as the real-world geometry remains unknown.Rather, our intention is to use these approaches to assess the variability of these fault maps and the faults they containthe better to assess the uncertainty in the array of fault maps that are produced.

Identifying, tracing, and mapping faults in seismic images
In this paper, we use the term "mapping" explicitly to mean the act of creating a three-dimensional interpretation of an individual fault or arrays of faults within the seismic volume.Horizontal visualisations of the 3D volume are called plan views (as opposed to vertical visualisations, which are termed profiles).At individual locations, faults are identified in these 2D visualisations (plan views or vertical profiles) as points (nodes) that, in these 2D visualisations, are linked as fault traces (sinuous or polylines embedded in 3D space).These fault traces in serial profiles or plan views are linked to create a mapped fault surface.
Fault identification, the interpretation of fault traces and the creation of 3D maps can be very challenging.Structures, such as channel margins, unconformities, igneous and salt bodies, also disrupt stratal reflections and, therefore, can be mistaken for faults (Baudon & Cartwright, 2008a, 2008b, 2008c;Delogkos et al., 2020;Faleide et al., 2021;Walsh et al., 1996).Therefore, it is important to define explicitly the criteria by which faults are identified.Faults rarely manifest as continuous reflections themselves; instead, seismic datasets often image faults as disturbances in the continuity of stratal reflectors (Dimmen et al., 2023;Gibson et al., 2005;Iacopini et al., 2016a;Misra and Mukheerjee, 2018).Analysing the continuity of stratal reflections and their terminations are the foundations of almost all fault-identification strategies.How faults manifest depends on two groups of factors: 1) the frequency content, the quality, the acquisition, and the processing characteristics of the data, and 2) the interpreter's criteria to recognise potential faults (Townsend et al., 1998).Seismic resolution exerts considerable constraint in fault recognition; faults below the seismic resolution cannot be detectable, or their manifestation is more challenging to identify (Dimmen et al., 2023).Lastly, the criteria employed can vary according to the interpreter's expertise and knowledge and the geological complexity of where the interpretation is undertaken.
Additionally, seismic attributes are excellent tools for highlighting discontinuities within seismic volumes (e.g., Kim et al., 2021).They enhance the variation of a particular seismic characteristic, such as amplitude, dip, azimuth, and curvature, that anomalies can be associated with the location of faults (e.g., Chopra and Marfurt, 2007;Libak et al., 2017).Their application generates a new volume visualisation (including vertical and plan views) that, most of the time, supports the interpretation of faults.Attribute anomalies are not necessarily fault expressions; they could be sedimentary structures or related to the presence of pore fluids.Their use requires a good understanding of attribute properties and care in adjusting the appropriate display setting (Marfurt and Alves, 2015;Townsend et al., 1998).Many seismic attributes are available and have proven to be practical tools for visualising and identifying faults in seismic volumes (e.g., Cohen et al., 2006;Kumar and Sain, 2018;Libak et al., 2017;Kim et al., 2021).Seismic attributes visualised on time slices are commonly used in association with vertical profiles to create fault interpretationsbut examples of fault-picking methodologies using these types of plan views alone are rare.
Vertical profiles are by far the most common ways of illustrating fault interpretations, and mapping tools in software packages are generally designed to use these views in building fault maps.The most distinct  criteria correspond to offsets in reflections, amplitude changes, abrupt changes in the dip of stratal reflections together with narrow monoclines (Iacopini et al., 2016a;Solum et al., 2016;Townsend et al., 1998;Misra and Mukheerjee, 2018, Fig. 2a).The accurate correlation of reflections between both sides of a potential fault is crucial because miscorrelation can drive errors in further analysis such as fault displacement profiles.In addition, display settings play a key role when using seismic profiles.For example, the employment of vertical exaggeration affects the accuracy of fault mapping, amplifying the issue of apparent continuity (Alcalde et al., 2019).By assuming a simple seismic velocity structure in our study, we have made structural interpretations on vertical profiles with minimal vertical exaggeration, an approach we recommend for all fault mapping on seismic images.
Moreover, interpreters often use time slices populated with seismic attributes (variance, dip) to analyse fault patterns in plan view (Fig. 2b).Also, they employ these views to correlate fault traces when using   (Bahorich et al., 1995;Pepper and Bejarano, 2005).
Finally, horizon maps are visual representations of stratal surfacesa basic product of a standard 3D seismic interpretation.Normal faults with significant heaves are represented by gaps in horizons.Smaller faults can be inferred from abrupt dip changes and narrow monoclines in mapped stratal surfaces (Fig. 2c; Dimmen et al., 2023).Alcalde et al. (2019) note that these approaches are sensitive to the vertical scaling of the seismic imageand vertical exaggeration can yield misleading results.Regardless of this effect, the use of horizon maps to identify, trace and map faults carries additional uncertainty linked to the interpretation and correlation of the stratal reflectors themselves.

Seismic image set, visualisation, and analytical platforms
Our study utilises an anonymous high-resolution 3D time-migrated seismic volume with low noise content.The seismic data were acquired with a bin size of 12.5 × 12.5 m, which is the spacing between Inlines and Crosslines.The data follow approximately zero-phase with the seafloor reflection displayed with negative polarity.The polarity colour scheme varies between blue and red, indicating a decrease and increase in impedance, respectively.Seismic volume frequency content ranges between 8 and 70 Hz (Fig. 3, b).We focused our interpretation on a shallow seismic package as uncertainty in fault placement increases with depth (Schaaf and Bond, 2019).Also, this seismic unit does not present a complex structural history, seemingly only displaying faults with normal offsets.The reflective sediment package comprises strong acoustical reflections with excellent lateral continuity without significant unconformities or complex depositional geometries, making the seismic succession ideal for structural studies (Fig. 3a).
Our work employed four software packages to access a range of analysis and interpretation tools.The interpretation process was carried out using Petrel E&P software (Schlumberger), measurements in MOVE Suite (supplied under academic licence by Petroleum Experts), the analysis of statistics in MINITAB (version 19) and the generation of fault maps in ESRI® ArcGIS Pro.

Fault interpretation workflows
We generated five fault maps by applying distinct workflows to evaluate the influence of methods in interpreting 3D seismic volumes.Fig. 4 contains sketches illustrating how each fault mapping was carried out, and the workflow steps are briefly described.Fig. 5 shows a table summarising the main approaches of each workflow.Our methods can be grouped into two categories based on the type of view where faults were mapped: i) plan view, and ii) vertical profiles.Workflows 1, 2 and 3 are part of the group i).Workflow 1 used plan view time slices draped with variance attributes.In Workflow 2, maps of stratal horizons were created by using vertical profiles to inform potential fault locations projected into plan views.Workflow 3 employed the same stratal horizons interpretations as Workflow 2, from these, surfaces are created and draped with the seismic attribute of variance, before fault interpretation.Workflows 4 and 5 are part of the second group (ii) and utilised vertical profiles, with differing sample densities.In this regard, Workflows 1 and 4 did not employ under-sampling, rather using every available image, be that every time slice (Workflow 1) or vertical profile (Workflow 4) -a time-consuming interpretational choice.In contrast, Workflows 2, 3 and 5 were under-sampled to replicate more time-efficient approaches more commonly adopted by interpreters.How fault placement is derived also varies according to these various methods.For Workflows 1 and 3, faults are placed along elongated shapes displaying high variance values.In Workflow 2, faults are defined by heave maps, the gaps left by horizon mapping.In Workflows 4 and 5, fault placement is determined by the recognition of reflection discontinuities.As seismic reflection data have a particular resolution, there is an inherent uncertainty in the placement of nodes and lines when interpreting faults.In this context, our measurements are limited by data resolution and precision, which varies between the workflows, this is discussed in more detail below.

Workflow 1. fault interpretation using plan view time-slices in conjunction with variance attribute
Seismic attributes, such as amplitude contrast enhancement, dip illumination or edge enhancement, are routinely employed to ease fault recognition and speed up fault interpretation.Their use needs to go hand in hand with a good understanding of their properties and awareness of potential pitfalls (Townsend et al., 1998;Marfurt and Alves, 2015).Variance attribute is widely used as a discontinuity-processing algorithm which efficiently enhances discontinuities in the continuity of reflections (Van Bemmel et al., 2000, Pepper andBejarano, 2005).These highlighted discontinuities, representing zones of low coherence, can represent faults (Gibson et al., 2005).In this workflow, we utilised every time-slice through the variance volume.It does not involve under-sampling.In each time-slice we interpreted faults in places where variance values were high, linking continuous traces as poly-lines.These poly-lines were then connected, moving through the stack of time-slices to create arrays of fault surfaces, using vertical sections as guides (Fig. 4).In this, we assume a positional uncertainty of at least 8 m, which is the minimum vertical resolution according to the analysed reflective package.

Workflow 2. fault interpretation in plan view using horizons picked on vertical profiles
Workflow 2 utilises the mapping of five stratal horizons in vertical seismic profiles to define fault gaps visualised in plan view.Five seismic reflections with high amplitude, including the seabed, were mapped across the full seismic volume using seismic profiles with amplitude display, and a gap was left where fault breaks were interpreted.Each of the five horizons was picked manually using every in-line, so with a sampling density of 12.5 m.Faults were identified where reflectors were broken with offset, abruptly deflected, or dimmed in amplitude.The result was five horizon maps with gaps representing fault heaves.Faults were created by placing nodes along the centre of the horizon-fault gap in plan view.The faults derived from this mapping comprise five layers of polylines.Where the polylines for different horizons clearly stack above each other, they can be interpolated to create a single fault  surface.In more complex areas, to avoid aliasing, the individual fault traces in the spaced horizons were tied using an in-line vertical profile (Fig. 4).The precision of the location of nodes is determined by the halfwidth of the fault gap, which is, on average, 15 m.

Workflow 3. fault interpretation in plan view using horizon surfaces created from the horizons interpreted in workflow 2 draped with variance attribute
Workflow 3 uses variance attributes draped onto horizon surfaces created from the horizons interpreted in Workflow 2. As in Workflow 1, nodes and lines were placed where high variance values formed an elongated shape representing zones of low coherence (high variance).Five layers composed of irregular polylines representing fault traces were generated.Where the polylines for different horizons clearly stack above each other, they can be interpolated to create a single fault surface.In more complex areas, to avoid aliasing, the individual fault traces in the spaced horizons were tied using an in-line vertical profile (Fig. 4).We assume an uncertainty of 15 m as in Workflow 2.

Workflows 4 and 5. fault interpretation using vertical seismic profiles
These workflows interpret faults on vertical profiles.These two workflows use the same basic mapping method.For a given vertical profile, fault identification involves the detection of disturbances in seismic reflection continuity.Fault mapping initially focuses on recognising reflection offsets, and nodes were placed in the middle of these offsets.Subsequently, the mapping is extended upward or downward by identifying other disturbances, such as reflection deflections or amplitude dimming.To correlate two fault traces found in two subsequent profiles, we required them to share a similar location, have a consistent character in the seismic image or disturb the same seismic reflections and have the same dip direction.Faults are mapped one by one through the seismic volume.When it becomes no longer possible to recognise and correlate a disturbance, the fault trace is ended (effectively identifying part of the fault tip-line), and the process starts again with another fault (Fig. 4).In Workflow 4, we carried out the fault mapping using the full image set (without under-sampling), and so has a profile spacing of 12.5 m.Workflow 5 utilised this full interpretation but was undersampled at a spacing interval of 50 m (one in four in-lines).Workflow 4 inherited the interpreted fault sticks from Workflow 5, every fourth line, and correlated between them.This approach was adopted to avoid introducing further uncertainty that can arise from creating multiple interpretations of the original seismic imagery.For Workflow 4, we consider the spatial resolution to be 6.5 m, half the spacing between Inlines.In the case of Workflow 5, we estimate this uncertainty to be 25 m which corresponds to half of the applied under-sampling.

Results: qualitative analysis of fault maps
The five different workflows generated distinctly different 3D fault maps, that we now compare.Beyond providing an illustration of the uncertainties inherent in the interpretation of even high-resolution, high-quality seismic imagery, the following results illustrate the importance of documenting interpretation workflows.There are (caption on next column) different ways of making comparisons of differences and similarities in our maps that we now consider.

Visualising on a common plan view
The fault maps created by the five workflows are compared in Fig. 6, viewed from above.The five maps display faults with a preferred orientation, striking NW-SE.Workflow 1 produced a total of 141 fault surfaces characterised by a dense number of polylines (Fig. 6a).The maps derived from Workflows 2 and 3 are formed by 73 and 43 fault surfaces, respectively (Fig. 6b and c).Workflows 4 and 5, based on profiles, generated a total of 82 and 75 fault surfaces, respectively (Fig. 6d and e).
We use two different types of display to illustrate qualitatively the differences in interpretation outcomes from the five Workflows.To analyse the fault pattern in two dimensions derived from our five maps, we chose the time slice with the highest number of faults interpreted in Workflow 1 and took it as a referential horizontal view.We then combined the other four interpretations at this referential horizontal view to obtaining their fault pattern at that level.Each interpretation provided a 2D fault map (Fig. 7).The visual examination of these fault maps reveals a marked similarity in fault patterns between workflows 2, 4 and 5.The map extracted from Workflow 1 stands out because it has the largest number of faults (100), and its fault pattern is characterised by more segmented small faults.On the other hand, Workflow 3 shows more continuous fault trends and is the map with the lowest number of faults (39).Despite these divergences, the five maps distinctly display the same preferred fault orientation (NW-SE), and faults are located at roughly identical locations.These facts validate the five maps in the sense that they are coherent and map broadly the same fault pattern.

Heat map
Here, we present a heat map (Fig. 8) as a representation of the variation in fault interpretations generated by the five workflows.Heat maps are excellent and widely used tools for visualising spatially distributed data, displaying the overall shape and concentration trends (e.g., Netek et al., 2018), and for collating multiple interpretations of image sets (e.g., Torvela and Bond, 2011).The heat map is a helpful visualisation tool to evaluate consistency between our 2D mapping workflows, to show the variability in fault interpretations and could be used to locate parts of seismic volumes that have a higher probability of containing faults.
Using ArcGIS Pro, we constructed the fault heat map by applying the Kernel Density tool to our five 2D fault maps composed of lines.The Kernel function calculates the area of influence of each feature (lines) based on an entered ratio (25 m) and then estimates density values.This ratio is based on the maximum estimated spatial uncertainty, which arises from Workflow 5. A key choice in creating heat maps is the colour range, hue, and intensity, which enhances the perception of map consistency and indicates different clustering densities in the data between maps.Heat maps emphasise hot spots as places of high clustering.Here we use a discrete colour scale, as we are only comparing 5 maps, with each colour representing a different degree of density between lines (faults).The heat map (Fig. 8) displays elongated shapes with colour variation, showing greater consistency (red-brown colour) and greater inconsistency (purple colour) between the five fault maps in a single  8 shows and reinforces the fact that our five 2D fault maps have high consistency and share similarities in terms of fault location and trend.

Juxtaposition map
The workflows generated five 3D fault maps that show coherence in fault trend and distribution when viewed in 2D plan views.However, there are differences in linkage, length, and number of faults.To perform a qualitative analysis of these variations, we present a juxtaposition map (Fig. 9) showing the explicit superposition between the five fault maps.We created a line with a ratio of 25 m around each fault in ArcGIS software.We again applied this ratio because it corresponds to the maximum estimated spatial uncertainty (Workflow 5).In this way, our five 2D fault maps are now composed of closed polylines.These five maps were then juxtaposed and visualised as a single composite map, shown in Fig. 9.This juxtaposition map shows zones where closed polylines coincide and differ.Differences in fault pattern and linkage are evident (Fig. 9a-b).In Fig. 9a, Workflow 3 defined a single fault, whereas the other workflows produced two faults.The superposition of the five 2D fault maps also displays variations in fault continuity and lengththat we consider later.

Fault shapes
While the juxtaposition map (Fig. 9) reveals differences in fault linkage and continuity, the different workflows examined here generate distinctly different shapes for individual faults.We illustrate this by selecting a single fault that has been mapped in all five workflows (Fig. 10).The visualisations in Fig. 10 show the composite apparent tiplines for the mapped fault, defined by the limiting nodes on the constituent polylines used in the interpretations.Therefore, they show the extent and shape of the fault surface derived from each of the five workflows.There are several striking, though perhaps unsurprising, features that arise through the comparison of the five workflows.
First, we compare the maps for the two workflows that use the full image set in the 3D seismicworkflows 1 and 4 derived from serial plan view time-slices and serial in-line vertical profiles, respectively (Fig. 10a  and d).Not surprisingly, these show complex, highly serrated tip-line geometries.Maps generated using plan views of the seismic imagery tend to generate apparent fault tip-lines that are highly serrated (Fig, 10a).These serrations are oriented parallel to the trend of the fault.In contrast, those fault maps generated using profiles generate serrated apparent tip-lines where the serrations are aligned with the dipdirection (Fig. 10c).Note that the fault has a greater extent in the dip directly mapped on vertical profiles, whereas it has a greater lateral extent in the interpretation derived from plan views.
Second, the interpretations derived from under-sampled imagery can be compared with those derived from the interpretation of the full image set.Unsurprisingly, the under-sampled versions have significantly simpler tip-line patterns because they are defined by fewer nodes than those interpretations using the full imagery.Consequently, the fault surfaces mapped using under-sampled imagery approximate more closely to ideal planes.This observation confirms those made in a similar study by Michie et al. (2021), albeit they restricted their mapping to vertical profiles alone.

Quantitative analysis of fault parameters
The similarities and differences described above are qualitative, therefore, tentative as they are based on the comparison of different visualisations: fault maps, 2D plan views and fault shapes.In this section, our qualitative data (five 3D fault maps) are transformed into quantitative data to numerically examine the impact of the different workflows on fault outcomes and to express their differences and similarities in numbers.We extracted two parameters: the maximum fault length and fault surface area (Fig. 11a and 12a).
Fault length and area are dimension features that define the size of faults in 2D and 3D, respectively.Each fault mapping provided two sets of measurements which were then contrasted by using basic statistical measurements (i.e., median, mean, etc.) and different graphical forms (Figs. 11 and 12).The statistical analysis performed in this work is elementary but very useful for revealing differences between our five fault maps.These assessments of dimensional parameters (area and length) are important for the further quantification of displacement patterns both on individual normal and or linked normal fault segments,   together with statistical representations of relationships between fault size and displacement (e.g., Kim and Sanderson, 2005;Torabi and Berg, 2011).However, we do not develop these considerations here.
The statistical assessment was performed in MINITAB v19 software.We extracted two fault measurement sets from each fault mapping: (1) maximum fault length and (2) fault surface area.Each set was treated as a sample.We first analysed each fault size parameter separately and then integrated and compared both together.Descriptive statistics were used to obtain the central tendency and distribution values of each sample.The mean and median define the central tendency for each sample.The median is described as the numeric value separating the higher half of a sample from the lower half.The mean corresponds to the arithmetic average.The interquartile ranges, skewness, and kurtosis are distribution measures that describe data distributions.Interquartile ranges divide the sample into four equal parts, meaning interquartile 1 contains 25% of the sample and so forth.Skewness measures the symmetry of a sample distribution based on a normal distribution curve (bell-shaped).Finally, kurtosis analyses the outliers of a sample which define the distribution pattern.This value is also obtained by comparing the sample distribution shape curve with a normal distribution shape curve (bell-shaped), defining where values are concentrated in the tails or around the mean.All these central tendency and distribution measures are tabulated in tables (Fig. 11b and 12b).
Additionally, we plotted measurement sets and their statistics in graphs with the aim of clarifying and integrating more effectively the dissimilarities and similarities between the five mappings (Fig. 11c-e and 12c-e).We used box plots to display the elementary statistical measurements (Fig. 11c and 12c).In this chart, the parameters described above (mean, median, interquartile ranges, outliers, and sample values) are visualised.Boxes indicate where 50% of the sample lies, corresponding to interquartile 1 and 3.The T-shaped whiskers represent the location of the maximum and minimum values excluding outliers (diamonds), which are measures positioned 1.5 times further away than the interquartile distance.The solid and segmented horizontal lines point out the median and mean values, respectively.The points displayed on the left of each box are the measures of each sample, and in darker colours are the corresponding outliers.
Samples distributions were displayed and compared in scatterplots (Fig. 11d and 12d).The X-axis shows all values in ascending order, and the Y-axis displays the number of faults.Finally, we constructed bar charts displaying the total sum for each parameter (Fig.s.11e and 12e).

Maximum fault length
The first parameter corresponds to the maximum horizontal length, which is a two-dimensional fault size measurement.Fault interpretations in this study are defined in 3D as a constellation of points which, in turn, are assembled in 2D fault traces, sinuous polylines (i.e., fault sticks).To measure the maximum horizontal length of a fault, we first projected all points into a horizontal plane.Then, we constructed a line based on the best Polynomial fit.The length of this best-fit line is what we refer as the maximum fault length (Fig. 11a).In this section, we used the term fault length to mean the maximum fault length.
We now compare the distributions and patterns of fault lengths derived from the five workflows.While Workflows 2, 4 and 5 exhibit roughly identical features in the three plots shown in Fig. 11, Workflows 1 and 3 display differences.However, in general terms, the five samples exhibit similar statistical attributes and distribution patterns with measures highly concentrated at low values.The scatterplot (Fig. 11d) clearly exhibits this similarity, the five curves show and start with steep slopes because roughly 80% of each measurement set falls within a very tight range of length values.Similarly, the five curves flatten when reaching the extreme values (outliers) represented by diamonds.The curves of Workflows 1 and 3 differ in their positioning because of their number of faults.Likewise, Workflows 2, 4, and 5 are in close proximity due to their similarity in the number of faults.This distribution pattern is also validated by positive skewness and kurtosis values, meaning that fault length measures have asymmetrical distribution patterns with high clustering at low length values and extremes (outliers) located at high values.
As outlined above, Workflow 1 and Workflow 3 stand out, displaying differences when comparing them with the other maps.Workflows 1 and 3 are based on the use of plan views populated with variance attributes, Workflow 1 used time-slices, whilst Workflow 3 surfaces.The fault mapping derived from Workflow 1 has the largest number of faults with short fault lengths.The sum of these lengths of short faults from Workflow 1 is almost twice that of the total summed lengths from the other four Workflows (Fig. 11e).It is worth remembering that Workflow 1 used the full image set.The fault map derived from Workflow 3 contains the longest fault length measurement, the greatest range, and the smallest number of faults interpreted.The total sum of fault lengths derived from Workflow 3 is roughly identical to the total sum of Workflows 2, 4 and 5, which contain twice the number of faults than the mapping derived from this method.
The numerical analysis corroborates the visual comparison of the 2D fault maps in Figs. 8 and 9: Fault maps derived from Workflows 2, 4 and 5 evidence close similarity.In contrast, fault maps derived from Workflows 1 and 3 moderately differ.Regarding Workflow 1: its fault map is characterised by the largest number of highly segmented faults with small fault lengths.Lastly, the fault map derived from Workflow 3 is distinguished by the smallest number of faults with more continuous and longer lengths.
Finally, the contrasting of all diagrams and tables shown in Fig. 11 permits us to establish that Workflows employing plan views (Workflows 1, 2 and 3) are more effective at interpreting faults with larger fault lengths than the methods using vertical profiles.

Fault surface area
The second parameter analysed is fault surface area, which is a threedimensional fault size measurement.All workflows define fault surfaces from the nodes on the polylines traced from the various workflows.For an individual fault, these nodes are points on the curvi-planar fault surface.However, various mathematic approaches exist for interpolating between these nodes to create the fault surface.These approaches, therefore, represent further diversity, beyond that represented by the (caption on next column) different mapping workflows developed here, in the models constructed for each fault.We reserve discussion of these interpolations for a future paper.Rather, we use part of the fault node population to illustrate quantitatively the differences between maps created from the five workflows.
We used the endpoints of each fault stick/polyline to create a continuous loop delineating the tip-line of a fault.These fault tip polylines were converted into polygons to obtain a fault surface.Polygons preserve the border of their preceding polylines and are curvi-planar surfaces (segmented line; Fig. 12a).We finally measured the area of these polygons to obtain the fault surface area.The analysis of these surfaces is a simplification, as we did not consider all middle nodes, which add roughness to surfaces.
Additionally, a single depth conversion was completed to extract fault surface area values, as the seismic volume is displayed in two-waytime.We assumed a constant velocity and performed a depth conversion in which the converted metric z-axis equals the time (1 msec = 1 m).As our choice of seismic velocity is arbitrary, the individual numbers are not significant, but as we use the same depth conversion for all outputs, the comparison of the fault surface area between maps is valid.
We now compare the distributions and patterns of fault surface areas derived from the five workflows (Fig. 12).The five mappings have roughly similar statistical attributes and distribution patterns, with fault area measures largely concentrated at low values.For example, the segmented and continuous lines representing mean and median parameters nearly overlap, positioned within a narrow range.The similarity in distribution patterns is clearly visualised in the scatterplot of Fig. 12d, where the five curves show and start with steep slopes because roughly 80% of each measurement set falls within a very tight range of area values.Similarly, the five curves flatten when reaching the extreme values (outliers) represented by diamonds.The curves of Workflows 1 and 3 differ in their positioning because of their number of faults.Likewise, Workflows 2, 4, and 5 nearly overlap due to their similarity in the number of faults.This distribution pattern is also validated by positive skewness and kurtosis values, meaning that fault area measures have asymmetrical distribution patterns with high clustering at low area values and extremes (outliers) are located at high values.
The sum bar chart (Fig. 12e) unveils a relationship between the number of faults and the total area mapped.The fault map derived from Workflow 1 has the highest number of faults, therefore, the highest total fault area mapped, while the fault map with the lowest number of faults (Workflow 3) captured the lowest total fault area.Moreover, the total fault surface area mapped derived from the two methods without undersampling (Workflows 1 and 4) reached similar values.
In particular, fault area measures derived from the mapping of Workflow 1 exhibit the most disproportional distribution pattern, with the largest number of faults clustered at low values and including the lowest minimum and highest maximum fault area values within the five samples.Skewness and kurtosis measurements of Workflow 1 (Fig. 12b) support this disproportionality as they are by far the highest.
Outlier values represented by diamonds in the box plot and scatterplot (Fig. 12c and d) also show some similarities.The highest extremes of the two workflows based on the interpretation plan views populated with variance attributes reach almost identical values (Workflows 1 and 3; purple and light blue).Similarly, the outliers of the two workflows based on profiles are similar (Workflow 4 and 5, red and green).The highest outlier belonging to Workflow 2 is located halfway between these two pairs (Workflows 1-3 and Workflows 4-5).It is worth remembering that Workflow 2 utilised profiles with amplitude display to generate stratal horizons (plan views) where faults were interpreted.

Comparing fault dimension features
We examined the impact of the five distinct workflows through the comparison of two features that define fault size: maximum fault length and fault surface area.We now integrate and explore relationships between these two numeric variables (length and area).In terms of similarities, the five fault mappings exhibit broadly similar statistical attributes (i.e., mean median and interquartile ranges) and their correlation is well visualised in the boxplot diagrams in Fig. 11c and 12c.Moreover, they consistently produce the same distribution patterns characterised by measures largely concentrated at low values with the presence of outliers at high values.The similitude in distribution patterns can be seen in the scatterplots in Fig. 11d and 12d.These two major similarities endorse the coherence between the five fault maps observed in the qualitative analysis of fault maps in section 5.1.
However, the five fault mappings are not identical in terms of fault length and fault area, and the comparison of these two features outlined differences.Mappings derived from plan views (Workflows 1, 2 and 3) show moderately more asymmetrical and scattered distributions than mappings derived from the interpretation of profiles (Workflows 4 and 5).On the other hand, outlier values, which are measures that significantly differ from the others, reach different values.Workflows employing plan views (Workflows 1, 2 and 3) contain extreme values with larger fault dimensions compared with the outliers of Workflows based on profiles (Workflows 4 and 5).Hence, the interpretation of faults using plan views is more effective in capturing larger fault sizes, in terms of fault length and fault surface area.
The two methods that did not involve under-sampling in their workflows have roughly identical total sums of fault surface area values but differ in their total sums of fault length values.This discrepancy can be related to the type of section employed; Workflow 1 utilises plan views (time slices) which optimise the mapping of horizontally longer fault traces.In contrast, Workflow 4 used vertical profiles undermining horizontal dimensions.In contrast, the two fault interpretation methods based on the interpretation of vertical profiles do not show considerable differences between their mappings, as Workflow 5 is the under-sampled version of Workflow 4. So, their comparison unveils the impact of undersampling.The interpretation of all vertical profiles captured a wider range of fault sizes, while the under-sampling restricted the range captured.Fault sizes derived from Workflow 4 are characterised by shorter and longer lengths and smaller and larger areas.Thus, the application of under-sampling reduces the size of the potential faults.
Finally, we emphasise the importance of the documentation of the methods, which allows us to explore and integrate the differences and similitudes between the five fault maps.Failing to document the mapping workflows would have prohibited such analysis, and therefore, any statistics derived from these maps have unquantifiable interpretation uncertainties.

Fault shape comparison
This section compares quantitatively the five versions of the fault shown in Fig. 10.This more detailed examination of a single fault considers the reasons behind the divergences or convergences between the five methods compared in this study.We aim to compare five measurements, incorporating the area and length measurements that have already been explained above.The three new measurements we now consider are the number of nodes, perimeter, and height.The number of nodes refers to the number of limiting nodes that compose the tip-line for each mapped fault.The perimeter is the length of the tip-line, and the height is the vertical difference between the uppermost and lowest nodes.All these measurements are reported in Fig. 13, and the five visualisations are in Fig. 10.
The sampling density determines the number of nodes of a fault.A fault will have a high number of nodes when the interpretation is carried out in a high number of sections.In decreasing order, Workflow 1 is the most densely sampled fault, followed by Workflow 4, then Workflow 2, Workflow 3, and finally Workflow 5 (Fig. 13; number of points).Unsurprisingly, the two workflows (1 and 4) that did not involve undersampling generated faults with the greatest number of nodes.
Fault perimeter, the length of the mapped fault tip-line, quantifies node alignment for a given fault area.The length of a polyline (perimeter) composed of node points aligned in a single plane will be less than one composed of spatially distributed node points.The biggest perimeter is produced by Workflow 1, followed by Workflow 4, Workflow 2, Workflow 5, and Workflow 3 (Fig. 13; perimeter).Our data evidence a direct relation between the number of nodes and the perimeter.Particularly, Workflow 3 shows a high alignment between points because even having more nodes than Workflow 5, it produces a smaller fault perimeter, meaning its 17 points are more aligned in 3D space than those fault nodes derived from Workflow 5.This is also consistent when looking at fault area values, the total area derived from Workflow 3 is roughly half of the total area produced by Workflow 5 (Fig. 13; area).
Regarding area values, the biggest area was produced by the application of Workflow 5, followed by Workflow 1, Workflow 2, Workflow 4, and Workflow 3. Area measures do not follow the same relative order as the number of points and perimeter because of the simplification process to create fault surfaces and the degree of alignment between points.
The contrasting of fault length measures in Fig. 13 endorses the remark noted in section 6: the utilisation of plan views generates interpretations of faults with larger lengths.By contrast, the comparison of height measures in Fig. 13 indicates that the utilisation of vertical profiles promotes the interpretation of faults with larger vertical extent.The last is visually evident in Fig. 10.Particularly, fault height derived from Workflow 3 is noticeably smaller compared with the other four measures.This anomaly is probably produced by the low sampling density of Workflow 3 and its higher susceptibility to drive misinterpretation.Lastly, the two methods that employ vertical profiles (Workflows 4 and 5) show similar fault heights, perhaps unsurprisingly, as both are derived from the same mapping.The number of nodes refers to the amount of ending points (tip points).Perimeter is the length of the polylines derived from joining all ending points using straight lines.The area is the area surface of the polygon embracing ending points.Length is the maximum horizontal obtained from the construction of the best-fit line based on the projection of all points onto a planar surface.Height is the difference between the uppermost and lowest points.

The impact of sampling density and type of sections in fault maps
The findings presented in the previous sections lead up to the following synthesis.Sampling density, the spacing between sections, and the type of sections, plan views or vertical profiles, are determinant choices in interpreting faults using seismic imagery.They significantly impact the dimensions and patterns of the resulting fault mappings.Fig. 14 summarises our results and shows a composite chart displaying the total sums for each workflow according to the following measurements: fault length, fault area and number of nodes.This figure integrates these three measures in a single plot (Fig. 14d), revealing the impact of following distinct workflows.
The five workflows analysed in this work captured fault mappings that roughly look similar in a quick visual examination.The heat map analysis evidenced consistency in fault trends.However, an examination in more detail revealed that the fault mappings have minor to significant differences in terms of fault length, fault area and fault patterns (Fig. 14d).Utilising every image available, either plan view or vertical profiles (Workflow 1 and 4), resulted in larger total sums of the number of nodes (Fig. 14a).Sparse interpretations led to smaller total sums of the number of nodes (Workflows 2, 3 and 5).Workflows employing plan views tended to capture faults with longer lengths than those based on vertical profiles.The involvement of under-sampling resulted in fault surfaces closer to perfect planes, therefore, smaller fault areas.On the other hand, methods employing the entire dataset, avoiding undersampling, resulted in more complex surfaces with larger fault area values.The fault patterns generated by workflows 1 and 3 stand out, evidencing distinct fault linkage.The fault map derived from Workflow 1 exhibits numerous short-length faults, and the one derived from Workflow 1 scattered long faults (Figs. 7 and 14d).In conclusion, the five workflows analysed in this study derived in coherent fault mappings, which in terms of fault length, fault area and fault patterns exhibit slight to considerable differences (Fig. 14).

Discussion
The five workflows employed to map the same 3D seismic volume generate different representations of the fault pattern in images.The consistency between these different outputs is demonstrated by the "heat map", in which the fault map outputs are stacked togetherrepresented in a single composite map view (Fig. 8).This visualisation is useful for showing, qualitatively, the parts of the seismic volume (or at least that part displayed in the map view) that are most likely to contain faults and the pattern of these most probable faults.This is important for grading areas in the subsurface for fault risk (e.g., for subsurface waste containment).However, the "heat map" obscures the variations in fault patterns and their dimensional parameters arising from the different workflows.The juxtaposition map emphasises mismatch, therefore, sites of greater structural uncertainty.Each map generates different expectations of fault continuity and linkage.But do these variations matter?Do different workflows create different fault patterns and geometries?And what are the implications of under-sampling seismic data when interpreting faults?The implications of the differences in fault pattern and geometry for interpretation workflow adoption and on further analyses and applications are discussed below.

The conventional use of vertical profiles
Although commercial pressures have driven the development of automatised workflows (e.g., Gibson et al., 2005;Kadlec et al., 2008;Pedersen et al., 2003), scientific publications are generally based on manual interpretations of 3D seismic volumes.Publications rarely include step-wise descriptions of how faults are mapped; consequently, assessing how fault outcomes are generated is unattainable.In our experience, many industry-based interpreters use a mixed workflow primarily founded on serial vertical profiles (our Workflow 5).Aliasing during the mapping can be minimised by consulting a selected time-slice displaying variance to track an individual fault along narrow bands of low seismic coherency.Additionally, interpretation based on vertical profiles is effectively promoted through the large number of figures that use this type of view to display mapped faults.Part of the reason for this may be the inheritance of 2D vertical profile displays from seismic imagery acquired in 2D, where interpretations relied entirely on vertical profiles.
We believe the community-wide reticence of using tiered plan views (rather than serial vertical profiles or a more balanced combination of the two) in seismic interpretation could limit the development of a better understanding of structures in the subsurface.Plan views are, of course, routinely used to map faults in landscapes and when using geodetic data (e.g., Carter et al., 2007;Delong et al., 2010;Donnellan et al., 2017;Johnson et al., 2014).Further, seafloor maps derived from high-resolution bathymetric data, including seismic imagery, are extensively utilised to map and display faults (Laor and Gvirtzman, 2023;León and Somoza, 2011;Somoza et al., 2021;Stuevold et al., 2003).

Implications in fault geometries, continuity, and linkage
Fault dimensional parameters are critical for a range of analyses in applied geology, such as seismic hazard analysis and the development of Fault Growth Models (e.g., Kim and Sanderson, 2005;Mark, 1977;Marret and Allmendinger, 1991;Walsh and Watterson, 1988).Fault dimensions determine fault connectivity vertically and horizontally.Displacement profiles are constructed to unveil fault propagation and slip history, in this sense, the geometry of fault surfaces is essential in studying the growth and evolution of isolated faults and fault arrays (e. g., Baudon and Cartwright, 2008a;2008b;2008c;Lohr et al., 2008;Long and Imber, 2010;Mansfield and Cartwright, 1996;Willemse, 1997;Ziesch et al., 2017).Additionally, there has been extensive research into the relationship between maximum displacement and fault length (e.g., Kim and Sanderson, 2005;Somerville, 2000).
The inherent uncertainty linked with this type of dataset and its interpretation, together with the fact that we do not know the real geometry of faults in the subsurface, put in doubt our understanding.Image resolution limitations play a crucial role in this sense because they limit and hamper the identification of faults in seismic images (e.g., Alcalde et al., 2017;Dimmen et al., 2023;Simm and Bacon, 2014).As well as the noise content generated by processing or depth-conversion methods makes the recognition of faults challenging.How faults manifest in seismic images depends on their scale, small-scale faults are more challenging because of their subtle manifestations (Dimmen et al., 2023).According to the method employed, small-scale faults can be missed and/or fault geometries can be misinterpreted, influencing individual faults or fault array geometries.The communication between faults is hampered as, according to the workflow, fault linkage patterns can vary between poorly to highly linked fault outcomes.The variation in fault patterns is evident in Fig. 7, where the 2D map derived from Workflow 1 shows a highly segmented fault array, while the fault linkage pattern obtained from Workflow 3 exhibits a more continuous and less segmented fault array.On the basis that we do not know the real geometry of faults in the subsurface, we cannot assume one fault network pattern is wrong, and the other is right.The application of multiple workflows offers the opportunity to evaluate different scenarios, thus handling the uncertainty better.
The five 2D fault patterns analysed in this work display broadly similar fault trends.However, they clearly differ in the number of faults interpreted and have differences in fault continuity.We can use the results of our analysis to discuss the benefits and limitations of these different workflows.What we found is that fault map interpretations based purely on plan view generated longer fault lengths, whilst vertical profiles resulted in faults with greater height.We propose that this results from the type of view utilised to carry out the interpretations.As the length is defined as the maximum horizontal distance, workflows employing plan views based on the generation of horizontal fault traces capture the maximum fault length unless under-sampled.In contrast, vertical profiles may mis-sample fault lengths due to the spacing of vertical profiles.Although, as height was measured in this work as the vertical difference between the uppermost and lowest nodes, methods using vertical profiles captured in more detail fault height variation and maximised fault height.Workflow 2, which combines horizon interpretations on vertical profiles with the interpretation of faults in plan view, shares the most similarities with the two workflows entirely based on the interpretation of vertical profiles, fault patterns and fault dimensional parameters.These similarities can be related to the use of vertical profiles to carry out the interpretation of horizons.This workflow falls in the middle of the fault statistics distribution.Michie et al. (2021), in their comparative interpretations of seismic profiles cut from a 3D volume, argue that better fault maps are created by under-sampled spaced vertical profiles than those using the full imagery.In our under-sampled workflows, we found that fault surfaces can hold less irregularities and rugosities because their construction is based on fewer nodes and lines.The significance of fault surface rugosity is unknown as we do not know the real geometry of faults in the subsurface.Thus, fault rugosity derived from seismic interpretation could be an interpretation artefact or real geology.The tendency to interpret simple fault geometries by under-sampling and relying on vertical profiles may not generate better outcomes, in other words, more accurate representations of faults in the subsurface.This may account for one of 'the five myths' concerning normal faults discussed by Ferrill et al. (2017a): normal faults are linear features in profile.Mechanical stratigraphy exerts significant control over the geometry of faults; dip variations or refraction occur in multilayer sequences because of competence variations (e.g., Peacock and Sanderson, 1992;Childs et al., 2009;Ferrill et al., 2017a,b) The efficacy of seismic attributes in fault interpretation has been widely analysed (e.g., Kim et al., 2021;Lou et al., 2019;Odoh et al., 2014;Qi et al., 2019;Zhang et al., 2014;Zheng et al., 2014).Marfurt and Alves (2015) brought to light potential pitfalls, emphasising that their employment can trigger misinterpretations depending on the interpreters' expertise.Another important factor when employing seismic attributes is the noise content.In this context, the 3D seismic dataset analysed in this work contains very low noise content, allowing the application of seismic attributes.In this study, two workflows included their use, particularly variance attribute (Workflows 1 and 3).The fault mappings derived from these two workflows exhibit general similarities with the other three based on amplitude displays.We hesitated to perform a comparison in more detail as it was not the scope of this work.However, this work proves that mapping faults in plan views draped with variance attribute generate coherent fault maps with the ones generated by interpreting faults using standard amplitude seismic displays.

Is it possible to define the best workflow?
Uncertainties in fault interpretation in seismic data, as highlighted here, make the determination of an ideal workflow intangible.But the value in understanding the range of fault networks and fault dimensional parameters derived from different workflows is considerable.A better understanding of how workflows impact fault outcomes (i.e., under-sampling or the type of views) is important for handling the uncertainty in fault interpretations and producing more accurate further structural studies.As shown in this work, distinct fault interpretation methods generated fault arrays with different fault patterns and geometries.We promote the implementation of multiple workflows, therefore, multiple fault outcomes as a strategy to handle better the uncertainty behind the interpretation of faults using 3D seismic volumes.
The discussion above emphasises that the interpretation of faults on seismic imagery carries inherent uncertainties.Whether these uncertainties are important depends on the purpose behind the mapping.For the traditional exploration of earth resources, seismic interpretation is likely to be only part of an investigation.Knowledge of the subsurface will be gained through the lifetime of the exploitation of the geological resources as further information is integrated with earlier seismic imagery.In this way, the structural interpretation can be tested, modified, or redrawn as further information comes to light.These incremental paths to subsurface knowledge are more problematic where they inform subsurface engineering projects directed at storage or disposal of waste materials.Given the general absence of workflow documentation in many scientific publications, it remains difficult to assess the uncertainties in our communal understanding of the geometry of fault structures derived from seismic interpretation.In the past, the interpretation of seismic volumes has been broadly restricted to deterministic interpretations due to the time and cost investments involved in these commercial engineering applications.Here we show how heat and juxtaposition maps created from different workflows might help to envision uncertainties in fault maps and inform risk associated with single deterministic fault models.
The inherent uncertainty involved in the generation of 3D seismic volumes, their data-rich nature, their large content of images, and the intrinsic variability of geology make each dataset unique.This uniqueness hampers the establishment of the 'best' workflow.The success of a determined workflow can vary depending on dataset particularities and uniqueness.This goes along with our intention not to recommend any of the workflows analysed here.Rather, our target is to enhance the variability of fault mappings generated by applying different workflows to the same image dataset.The qualitative and quantitative analyses of the five fault outcomes derived from the application of distinct workflows permitted us to establish three main workflows' particularities that are critical in determining the resulting fault outcome: (1) fault recognition and picking criteria, the principles that interpreters will use to track faults, (2) type of views and how they are blended, and (3) sampling density as it has been revealed by Cunningham et al. (2021); Michie et al. (2021); Ze and Alves (2016), 2019.

This work's limitations and future work
The apparent simple structural setting of the seismic volume analysed in this study plays an important role and limits our findings.Our five mappings showed that the 3D seismic volume analysed does not hold a complex structural history, evidencing only the development of normal faulting.This restricts our work to the tracking of normal faults.Even more, for small-scale normal faults, the faults interpreted in this work are not regional structures.Further studies are required for other and more complex structural settings, such as thrust belts or strike-slip faults.We emphasise that the abovementioned differences and similarities between the five fault mappings are particular to the image set analysed in this work.We aimed to show that different methods led to a range of distinct mappings.This approach is helpful to assess better the inherent uncertainty behind the interpretation of seismic imagery and compare multiple interpretation scenarios.

Fault interpretation and emerging challenges
Forecasts of the geometry of faults in the subsurface impact the assessments of other geological behaviours, such as cross-fault fluid transmissibility, as discussed by others (e.g., (e.g., Braathen et al., 2009;Childs et al., 2009;Marchal et al., 2003;Peacock, 2002;Walsh et al., 2003;Wibberley et al., 2008).Much of this activity, and therefore the common workflows in seismic interpretation, relate to the exploration and production of hydrocarbons.In recent years, new demands have emerged mainly in response to the challenges linked with industries of the Energy Transition.More detailed and accurate interpretations of the subsurface are required (e.g., Han et al., 2021).Advances and improvements in the interpretation of faults in seismic volumes are imperative, and methods are an essential starting point.Large datasets such as any subsurface seismic volume and the unknown real-world subsurface geometry require better and more sophisticated management of the uncertainty from the perspective of developing industries in geo-storage, for carbon, hydrogen and geothermal Energy, and the geological disposal of nuclear waste.The generation of multiple fault outcomes (stochastic models) by applying distinct methods offers many advantages to improve prediction and risk assessment by embracing uncertainty and calculating the likelihood of outcomes between several scenarios.

Conclusions
Almost all the knowledge about fault geometries in the subsurface derives from the interpretation of seismic images.This work has shown how different fault interpretation workflows, including different types of sections (plan views and vertical profiles), sampling densities and fault placement, lead to distinct 3D fault mappings.We compared qualitatively and quantitatively five fault mappings derived from five different workflows.Our findings are.
• 2D fault maps derived from each 3D mapping display broadly similar fault trends and distributions, indicating coherence in the way that the application of different workflows picked roughly the same faults.• We showed the efficacy of using a heat map to highlight the variability in fault interpretation generated by overlapping the five fault maps.It reveals interpretation consistencies/inconsistencies between the mappings and could be used as a likelihood map.• The employment of a juxtaposition map, showing the superposition of the five 2D fault maps, clearly uncovers differences in fault linkage and continuity, specifically in length.• Juxtaposition and heat maps as approaches for the integration of multiple interpretations and their comparison are excellent visualisation tools.• The comparison of a single fault evidenced that the different workflows applied here generate distinctly different shapes for individual faults.• Methods that use the full image set produce highly serrated tip-line geometries.This serration geometry varies according to the type of section employed.The interpretation of faults using plan views (time slices) generates faults with lateral serrated geometries, while interpretations derived from profiles produce faults holding vertical serrated geometries.• The application of under-sampling generates simpler tip-line patterns deriving in fault surfaces more closely to ideal planes.• The analysis of fault dimensional patterns, area and length, discloses differences associated with the particularities of each workflow.• Methods based on plan views produce fault surfaces with bigger areas and long lengths.• Sampling density is a critical factor in controlling the geometries of fault surfaces.It also determines the number of faults and total area mapped.More faults with greater fault dimensional attributes are interpreted when employing the full image set than when applying under-sampling.Workflows that involve high fault under-sampling may lead to fault misinterpretations.
Finally, we highlight the importance of understanding uncertainties in fault networks and geometries resulting from different interpretation workflows.Uncertainties in fault interpretation using seismic imagery,

Fig. 1 .
Fig. 1. 3D Seismic volumes. A. 3D visualisation of a Seismic Volume displaying standard seismic amplitude.B. 3D visualisation of a Seismic Volume displaying standard seismic amplitude.The segmented red line represents the interpretation of a potential fault in plan view and vertical profile.

Fig. 2 .
Fig. 2. What do faults look like? A. In seismic profiles, faults are identified through changes in the reflection's continuity, such as (1) variations in the reflection's dip (monoclines), (2) offsets and (3) amplitude dim.B. In plan view -time slices displaying coherence attributes such as variance.High variance anomalies can be interpreted as potential faults.C. In plan viewhorizon surfaces, faults are recognised by abrupt changes in the dip of horizons, these are then used to define faulty heave, seen as black gaps in the image.

Fig. 3 .
Fig. 3. Dataset. A. Seismic profile showing the lateral continuity of the reflective sediment package analysed.B. Data frequency content of the seismic volume utilised.The interpretation of faults is preferably carried out using imagery data set with high frequency.In this work, the reflective sediment package has a frequency of 50 Hz on average.

Fig. 4 .
Fig. 4. Workflow diagram.The Workflows are divided into two groups i) and ii).Group i) the interpretation of faults is based on plan views (dark grey background), and group ii) fault interpretations are based on vertical profiles (light grey background).Workflow 1: Fault interpretation is entirely based on the use of plan view time slices displaying variance attributes.Faults are placed at high variance values and later correlated.Workflow 2: Fault interpretation is carried out in plan views at horizon-fault gaps, using horizons interpreted on vertical profiles.This workflow uses vertical profiles at the early stage, and plan views are used to carry out the interpretation of faults.Workflow 3: The horizon interpretations from Workflow 2 are interpolated into horizon surfaces draped in variance attributes and visualised in plan view to map faults.Faults are placed at high values of variance.This workflow uses vertical profiles at the early stage, and plan views draped with variance attributes are used to carry out the interpretation of faults.Workflow 4 uses all vertical profiles available.Faults are tracked using seismic profiles and placed along reflection disturbances.Workflow 5 is an under-sampled version of Workflow 4. Faults are tracked using seismic profiles and placed along reflection disturbances.

Fig. 5 .
Fig. 5. Workflows features.Table showing the main characteristics of each workflow.

Fig. 6 .
Fig. 6.Plan view of 3D fault map derived from each Workflow.A. Mapping derived from Workflow 1: 3D mapping resulting from the interpretation of plan view time slices displaying variance attribute B. Mapping derived from Workflow 2: 3D mapping resulting from the interpretation of horizons in plan view containing fault gaps.C. Mapping derived from Workflow 3: 3D mapping resulting from the interpretation of plan view horizon surfaces displaying variance attribute.D. Mapping derived from Workflow 4: 3D mapping resulting from the interpretation of vertical seismic profiles.E. Mapping derived from Workflow 5: 3D mapping resulting from the under-sampling of the 3D mapping from Workflow 4. Group i): Workflows 1, 2 and 3. Group ii): Workflows 4 and 5.
Fig. 7. Visualisation of the faults in map view generated by each Workflow.The intersection of each 3D fault mapping with a horizontal time slice.This horizontal time slice corresponds to the time slice with the highest number of faults interpreted.A. 2D fault map view derived from Workflow 1. B. 2D fault map view derived from Workflow 2. C. 2D fault map view derived from Workflow 3. D. 2D fault map view derived from Workflow 4. E. 2D fault map derived from Workflow 5. Group i): Workflows 1, 2 and 3. Group ii): Workflows 4 and 5.

Fig. 8 .
Fig. 8. Heat map of the faults interpreted by all five Workflows.This map displays a composite view of all five fault maps generated by the different workflows.Fault density or overlapping degree of the five fault interpretations is represented as colours.Red-brown and purple represent high and low densities, respectively.

Fig. 9 .
Fig. 9. Juxtaposition map of the five Workflows.This map shows the superposition of the five plan view fault maps as elliptical polylines.Each polyline represents the location of a fault.Due to the high degree of overlap, we varied the visualisation of the different maps; we kept the polylines for the maps of Workflows 1, 2 and 5 and created coloured polygons for the maps derived from Workflows 3 and 4.

Fig. 10 .
Fig. 10.Visualisations of the interpretations of a single fault from each workflow.Each visualisation has in the background a seismic profile with a time slice displaying variance attributes.The boxes contain the polygon derived from joining the ending points of each interpretation, to create a visualisation of a single fault interpreted by the five methods.A. Fault visualisation resulting from Workflow 1. B. Fault visualisation resulting from Workflow 2. C. Fault visualisation resulting from Workflow 3. D. Fault visualisation resulting from Workflow 4. E. Fault visualisation resulting from Workflow 5.

Fig. 11 .
Fig. 11.Maximum fault length.A. Sketch showing how fault length values were obtained.B. Descriptive statistics.The table contains the values of mean, total sum, minimum value, Interquartile 1 (Q1), median, Interquartile 3 (Q3), maximum value, skewness and kurtosis.C. Fault Length Box Plot.This chart displays the central tendency and distribution patterns of each sample graphically.Coloured boxes represent half of each sample (50%).T-shaped whiskers indicate samples' ranges excluding outliers.Segmented and continuous lines represent mean and median parameters, respectively.D. Fault length scatterplot.The X-axis shows all values of fault length in ascending order, and the Yaxis displays the number of faults.E. Fault length bar chart.This plot displays the total sum of all length values for each workflow.

Fig. 12 .
Fig. 12. Fault surface area.A. Sketch showing how fault area values were obtained.B. Descriptive statistics.The table contains the values of mean, total sum, minimum value, Interquartile 1 (Q1), median, Interquartile 3 (Q3), maximum value, skewness and kurtosis.C. Fault Length Box Plot.This chart displays the central tendency and distribution patterns of each sample graphically.Coloured boxes represent half of each sample (50%).T-shaped whiskers indicate samples' ranges excluding outliers.Segmented and continuous lines represent mean and median parameters, respectively.D. Fault Area scatterplot.The X-axis shows all values of fault area in ascending order, and the Y-axis displays the number of faults.E. Fault Area bar chart.This plot displays the total sum of all area values for each workflow.

Fig. 13 .
Fig. 13.Comparison of fault parameters for a single fault.This chart displays the dimension parameters of the five different versions of the same fault.The number of nodes refers to the amount of ending points (tip points).Perimeter is the length of the polylines derived from joining all ending points using straight lines.The area is the area surface of the polygon embracing ending points.Length is the maximum horizontal obtained from the construction of the best-fit line based on the projection of all points onto a planar surface.Height is the difference between the uppermost and lowest points.

Fig. 14 .
Fig. 14.Summarising diagram.Composite chart integrating the total sum of fault length values, fault area values and the number of nodes of each workflow.A. Bar chart displaying the total sum of the number of nodes per workflow.Yellow starts indicate the positioning; for example, number 1 lies on Workflow 1 which is the method with the denset sampling method.B. Bar chart displaying the total sum of fault length values per workflow (identical to Fig. 11e).C. Bar chart displaying the total sum of fault area values per workflow (equivalent to Fig. 12e).D. Diagram showing the intersection between the total sum of fault length and the total sum of fault area.The yellow arrows indicate the direction in increasing order of the total sum of the number of nodes.The origin is Workflow 5 because it is the method with the most sparse sampling.This diagram also includes the fault pattern characteristics of the 'two extreme' fault mappings (Workflow 1 and Workflow 3).

Fig. 15 .
Fig. 15.The impact of methods on the interpretation of faults using seismic imagery.It illustrates the variability of fault interpretation controlled by the type of view employed.The multiple choices available to carry out fault interpretation, plus as we do not know the real geometry of faults in the subsurface, make this type of study highly uncertain.