How do differences in interpreting seismic images affect estimates of geological slip rates?

Uncertainties of geological structural geometry constructed based on seismic reflections can stem from data acquisition, processing, analysis, or interpretation. Uncertainties arising from structural interpretations and subsequent estimates of geological slip have been particularly less quantified and discussed. To illustrate the implications of interpretation uncertainties for seismic potential and structural evolution, I use an example of a shear fault-bend fold in the central Himalaya. I apply a simple solution from the kinematic model of shear fault-bend folding to resolve the geological input slip of given structure and then compare the result with a previous study to show how differences in structural interpretations could impact dependent conclusions. The findings show that only a little variance in interpretations owing to subjectivity or an unclear seismic image could yield geological slip rates differing by up to ∼ 10 mm yr⁻¹, resulting in significantly different scenarios of seismic potential. To reduce unavoidable subjectivity, this study also suggests that the epistemic uncertainty in raw data should be included in interpretations and conclusions.

1 Introduction

Geological studies typically integrate various spatial and temporal data that may have large uncertainties; these uncertainties originate from a lack of guiding models or accessible records in the natural settings, sometimes from necessary simplification in models, or from our limited understanding of how to interpret raw data. Also, given that geology is hermeneutic in nature (Frodeman, 1995), to reduce subjectivity, identifying, including, and quantifying uncertainties of interpreted information have been widely recognized as necessary (e.g. Wellmann and Caumon, 2018).

Two main foci of structural geology are evaluating reserves of natural resources and seismic hazards, which are respectively based primarily on established geological models as well as geological slip and its rate. Geological models are used to present spatial relations between faults and stratigraphic horizons; uncertainties of geological models have been widely studied and expressed (e.g. Lindsay et al., 2012; Thiele et al., 2016; Wellmann and Caumon, 2018; Randle et al., 2018; Schaaf and Bond, 2019). Geological or geomorphic slip, a metric of estimating seismic potential, is constrained by offsets or deformation in surface features (e.g. fluvial terraces, fold scarps) together with the proposed underlying structural geometry; their uncertainties have also been considered in several studies (e.g. Thompson et al., 2002; Davis et al., 2005; Bird, 2007; Zechar and Frankel, 2009; Amos et al., 2010; Li et al., 2013; Xu et al., 2020). Subsurface fault geometry – a critical element for constraining geological slip – is mostly simplified as a single planar fault, which can thus lead to large uncertainties (e.g. Brandes and Tanner, 2014; Butler et al., 2018; Drukpa et al., 2018). In some cases (when funding is less of a hurdle and field conditions allow it), the geometry is constrained by seismic reflection imagery, which is considered to yield a precise fault geometry. However, seismic reflection imagery may also suffer from several sources of uncertainties including the resolution of seismic images (e.g. Faleide et al., 2021), choice of velocity models (e.g. Totake et al., 2017), and bias in interpretations.

Figure 1 Two end-member models for shear fault-bend folding; after Suppe et al. (2004). (a) Simple shear and (b) pure shear fault-bend folding. The relations of the four geometric features – backlimb dip (δ_b), ramp dip (θ), backlimb length (L), and input displacement (d) – are described by Eq. (1).

Generally, a wide range of interpretations of a seismic image is possible mostly because seismic images are rarely clear enough to determine fault geometry conclusively. Moreover, differences in interpretations could be associated with an interpreter's prior knowledge and experience or guiding models with different underlying assumptions (e.g. Rankey and Mitchell, 2003; Bond et al., 2007, 2012, 2015; Groshong et al., 2012; Bond, 2015; Alcalde et al., 2017, 2019; Butler et al., 2018). The further impact of such differences in interpretations remains less discussed, especially how these differences could impact estimates of geological slip rates and thus seismic hazards and structural evolutions.

Here, I illustrate such an issue with a case of an unclear seismic image of a shear fault-bend fold in the central Himalaya. Firstly, I highlight a simple solution derived in Suppe et al. (2004). This solution can directly acquire the input displacement accommodated by shear layers in shear fault-bend folds based on the features in the shallow part of seismic images. Thus, the solution could be used to acquire total slip even when the quality of the seismic image is poor. This solution is validated with a case of a better-imaged structure – the Tainan anticline in the fault-and-thrust belt of Taiwan. Subsequently, I apply this solution to the case of the Main Frontal Thrust in the frontal part of the central Himalaya and compare the acquired total geological slip with a previous study. Finally, I quantitatively show how this difference in seismic image interpretation could impact the following conclusions regarding geological slip rates and thus implications for seismic potential and structural evolution.

2 Method

2.1 Geometry and kinematics of shear fault-bend folding

Shear fault-bend folds are commonly identified in thrust belts worldwide (e.g. Moore et al., 1990; Suppe et al., 2004; Shaw et al., 2005; Corredor et al., 2005; Yue et al., 2011; Almeida et al., 2018; Le Béon et al., 2019), and they are distinctly characterized by long, gentle, and rotated backlimbs. These characteristics are caused by shear strain in décollement layers (Suppe et al., 2004; Shaw et al., 2005; Hardy and Connors, 2006). Based on shear strain mechanisms, two end-members of kink-style kinematic models for shear fault-bend folding have been proposed (Suppe et al., 2004). First, the simple shear fault-bend fold has externally imposed bed-parallel and homogeneous simple shear accommodated by the basal shear layer, with no basal fault slip (Fig. 1a). Second, the pure shear fault-bend fold has no simple shear within the décollement layer (Fig. 1b); instead, its shortening is evenly accommodated by pure shear within the décollement layer, which moves along the basal fault.

Suppe et al. (2004) derived geometric relations among features in shear fault-bend folds of both types based on assumptions of conservation of area, line length, and bed thickness in the plane profile. Since these geometric relations are powerful for predicting unknown geometry based on partially revealed features, they are widely and successfully utilized in interpreting seismic reflections.

2.2 Acquiring total slip of given structures from seismic images

For tackling issues of seismic potentials and structural evolutions, the most critical element needed is the total slip along the targeted structures. Given that the input displacement along the décollement (i.e. input slip or shortening where décollement is horizontal) is mostly nearly horizontal, the total displacement cannot be identified directly from seismic reflections. Therefore, the slip amount is mostly indirectly decided based on applications of area–depth strain (ADS) analysis, or the area-of-relief method (e.g. Epard and Groshong, 1993, 1995; Gonzalez-Mieres and Suppe, 2006, 2011; Schlische et al., 2014; Groshong, 2015; Eichelberger et al., 2017; Wang et al., 2018). This method, assuming that area is conserved along the profile, requires at least two pairs of area–relief (i.e. excess areas of a horizon and decided reliefs from the given horizon). However, these features cannot always be determined in an ambiguous seismic image or could be incorrectly determined in a structurally complex area (Carboni et al., 2019), potentially leading to a large uncertainty in estimates of geological slip along a given structure.

Here, I apply geometric derivations from Suppe et al. (2004) and show that while this solution does not require the ramp locations or area–relief pairs in seismic images to be determined, it can independently constrain the slip along the décollement based on geometric features which are mostly clearer than other features in seismic images.

According to the geometric solutions in Suppe et al. (2004)'s shear fault-bend fold model, the input slip (d) accommodated by the shear strain within the décollement layer above the basal fault can be decided by Eq. (1) (Eq. 6 in Suppe et al., 2004):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\text{cos}\left(\frac{{\mathit{\delta }}_{\mathrm{b}}}{\mathrm{2}}-\mathit{\theta }\right)}{d}}={\displaystyle \frac{\sin \left(\mathit{\theta }\right)}{\mathrm{2}L\sin \left(\frac{{\mathit{\delta }}_{\mathrm{b}}}{\mathrm{2}}\right)}}.\end{array}$$

This equation determines the relations between the input slip (d), representing the total slip of the structure here, ramp dip (θ), backlimb dip (δ_b), and the length between the two axial surfaces that define the backlimb length (L).

Applying this solution to decide input slip has three advantages. First, L and δ_b specifically reflect the amounts of shear strain in the décollement layer (Fig. 1), and generally, due to their shallower positions, they could be more easily determined in the seismic images than other geometric features. Second, θ is a stratigraphy-controlled factor, so if it is unclear in seismic images, its possible range could be determined based on regional geological information. Third, the input slip (d) obtained by this solution is independent of types of shear strain (i.e. pure or simple shear), which are mostly indistinguishable in seismic images. With the deterministic relationship between these four variables (L, θ, δ_b, d), d can be assigned as the response variable and the other three variables (L, θ, δ_b) as the predictor variables. These predictor variables can be assigned reasonable values by referring to seismic images or regional geological information.

2.3 Including and propagating uncertainties

To minimize subjectivity in interpreting seismic images, I include uncertainties inherited from ambiguous seismic images and decide a range of possible realizations for each variable, which might be interpreted differently by individuals, together with assignments of the probability distribution functions (PDFs) for the predictor variables. These PDFs represent our best knowledge of these variables currently, and their widths of distributions could be narrowed down once the fault geometry is better resolved. Then, uncertainties from predictor variables propagate to the response variable by Monte Carlo simulations.

3 Applications

I firstly validate the highlighted simple solution by applying it to a case of the shear fault-bend fold in the foreland of Taiwan (Le Béon et al., 2019). Next, I apply this solution to an analogous case in the central Himalaya (Almeida et al., 2018), and then show the differential estimates of total slip and further implications originating from the difference in interpretations.

3.1 The Houchiali backthrust in the fold-and-thrust belt of Taiwan

The Tainan anticline has been categorized as a pure shear fault-bend fold above the Houchiali backthrust in the fold-and-thrust belt of southwestern Taiwan (Le Béon et al., 2019). Le Béon et al. derived total accommodated shortening according to the ADS analysis (e.g. Groshong, 2015); in their analysis, least squares linear regression of pairs of depths of interpreted horizons and corresponding excess areas are used to constrain the total shortening.

Figure 2 (a, b) The interpreted and uninterpreted seismic profile – D5, Fig. 2 in Le Béon et al. (2019). (c) The interpretation of this study and the decisions of predictor variables (i.e. L – backlimb length, θ – ramp dip, δ_b – backlimb dip). The pre-growth strata are interpreted in cyan, and the growth strata are interpreted in lavender. The distributions of predictor and response variables are shown in Fig. 3.

Figure 3 The distributions (a–c) of the predictor variables (L – backlimb length, θ – ramp dip, δ_b – backlimb dip) and (d) the response variable (d – total slip) for the Houchiali backthrust. The total slip estimated by this study (${\mathrm{0.522}}_{-\mathrm{0.049}}^{+\mathrm{0.050}}$ km) and Le Béon et al. (2019) (0.529 ± 0.044 km), shown by red fill, almost fully overlap.

3.1.1 Deciding the predictor variables

Independently from Le Béon et al. (2019)'s interpretation, I assign predictor variables of backlimb length (L), which is the distance between two axial surfaces (AS1 and AS2 in Fig. 2c) that bisect the interlimb angles in the theoretical model, as a uniform distribution with a range of ∼ 1.9 to 2.1 km. The ramp dip (θ) is assigned to be normally distributed with a mean of 29^∘ and a standard deviation of 1^∘. The backlimb dip (δ_b) of pre-growth strata (i.e. reflectors interpreted in cyan) falls within ∼ 7 to 9^∘ (Fig. 2c). I then apply the Monte Carlo simulation with 1 million random samples from the space of predictor variables, distributed according to the assigned PDFs, to simulate the full range of possible realizations of the predictor variables and arrive at the distribution of the response variable (d) (Fig. 3).

3.1.2 Results

The total slip estimated here (${\mathrm{0.522}}_{-\mathrm{0.049}}^{+\mathrm{0.050}}$ km) is nearly identical to Le Béon et al. (2019)'s result (0.529 ± 0.044 km). This similarity shows that nearly identical values can be acquired through different methods if the seismic image is properly interpreted. Also, the similarity demonstrates the effectiveness of the proposed solution to resolve slip along a given structure with a single planar ramp. The proposed solution is procedurally efficient, and not requiring certain steps, such as deciding area–relief pairs, reduces potential artificial uncertainties. Furthermore, two of three decided features could be defined from the shallow part of seismic reflections, and parameterized features can be processed without visualizing the full geometry.

3.2 The Bardibas thrust in the front of the central Himalaya

Next, the example of the Bardibas thrust in the foreland of the central Himalaya, Nepal, is examined. The Bardibas thrust is the foremost strand of the Main Frontal Thrust (MFT), which is suggested to accommodate most of the shortening across the Himalayan wedge in the Holocene (e.g. Lavé and Avouac, 2000), and thus the concluded slip rate of the Bardibas thrust has implications for the seismic potential in the Himalaya. The fold above the Bardibas thrust has been identified as a pure shear fault-bend fold, and the total slip accommodated has been decided based on only one pair of excess area and decided depth of décollement. The décollement depth was stated by Almeida et al. (2018) to have been based on Suppe et al. (2004)'s shear fault-bend fold model together with a decided ramp location. However, first, the quality of the image is not high enough to decide the mentioned features conclusively. Second, we cannot rule out the possibility of an existing upper décollement, along which mass can flux out; in such a case, the area used to derive the total slip, based on a calculation from a pair of area–relief, can be underestimated and the total slip is thus overestimated. Third, the shear mechanism for the shear fault-bend fold is decided to be a totally pure shear by Almeida et al. (2018), but if any simple shear involves the shear layer, the total slip – according to their calculation – would be underestimated (e.g. Epard and Groshong, 1993).

Here, independently from the decided features, I apply the simple solution proposed and validated in the previous sections to acquire the total slip of the Bardibas thrust. I then compare the result with Almeida et al. (2018)'s interpretation and discuss further implications.

Figure 4 (a, b) Interpreted and uninterpreted seismic profile; Fig. 2 in Almeida et al. (2018). (c) Assigning the ranges of possible values for predictor variables (i.e. L – backlimb length, θ – ramp dip, δ_b – backlimb dip) according to Almeida et al. (2018)'s seismic image of the Bardibas thrust. Their distributions are shown in Fig. 5.

3.2.1 Deciding the predictor variables

The anticline above the Bardibas thrust is in fact not exclusively a shear fault-bend fold (Suppe et al., 2004), but it is affected by fault-propagation folding (Suppe and Medwedeff, 1990). Fault-propagation folding specifically can complicate the length of the backlimb (L), which is used to indicate the shear strain in the basal shear layer. Where a fold is a combination of both shear fault-bend folding and fault-propagation folding, L may not straightforwardly indicate the strain in the shear layer. In this structure, the unconstrained portion, which is attributed to fault-propagation folding, of the L is treated as an epistemic uncertainty that is related to the quality of the seismic image. Thus, I include a possible range of L, inferring shear strain in the décollement layer. The minimum length (L_min) is 2.7 km, which is the length between axial surfaces 1 and 2 (AS1 and AS2), where bedding starts to change its trend. The maximum length (L_max) is 3.8 km, which is the length between axial surface 1 and 3 (AS1 and AS3), where the bedding starts to flex downward (Fig. 4). Thus, I assign L as uniformly distributed in a range of ∼ 2.7 to 3.8 km. The backlimb dip (δ_b) is not changed by the evolution of the fault-propagation folding and it is relatively clear; I assign it a range of 10–12^∘ with a uniform distribution. The ramp location and ramp dip (θ) could not be conclusively determined from Almeida et al. (2018)'s seismic image. Referring to regional bedding attitudes (e.g. Hirschmiller et al., 2014; Shrestha et al., 2019), I accordingly assign a normal distribution, with a mean of 30^∘ and a standard deviation of 5^∘ (Fig. 4), and assume it is a single planar ramp.

Figure 5 The distributions (a–c) of the predictor variables (L – backlimb length, θ – ramp dip, δ_b – backlimb dip) and (d) the response variable (d – total slip) for the Bardibas thrust. The total slip estimated in Almeida et al. (2018, 1.67 km), shown by the vertical red line, varies significantly by ∼ 0.6 km from the most probable value (${\mathrm{1.07}}_{-\mathrm{0.16}}^{+\mathrm{0.36}}$ km) estimated in this study.

Figure 6 Hypothetical ages of the Bardibas thrust (horizontal axis), the corresponding fractions of the shortening rates across the central Himalaya (right axis, overall shortening rate is ∼ 19.4 mm yr⁻¹ based on Stevens and Avouac, 2015), and the corresponding differential rates (left axis). These differential rates originate from the varying estimates of the total slip by Almeida et al. (2018) and this study. The black dots show the most probable values, and the grey shaded area shows the 68 % confidence interval.

3.2.2 Results and further implications

The total slip of the Bardibas thrust is estimated in this study to be ${\mathrm{1.07}}_{-\mathrm{0.16}}^{+\mathrm{0.36}}$ km, which is within the 68 % confidence interval; this is significantly lower than ∼ 1.7 km inferred by Almeida et al. (2018) (Fig. 5). I examine here how the different estimates of the total slip can impact the following derivations of geological slip rates, which serve as metrics to assess seismic potential when compared with the geodetic shortening rate across the given structure. Taking the estimated total slip ∼ 1.07 km and dividing it by the expected geological slip rate, the Bardibas thrust should be activated at ∼ 61 ka. The expected geological slip rate is suggested according to the geodetic shortening rate across the central Himalaya (∼ 19.4 mm yr⁻¹, Stevens and Avouac, 2015). If the thrust has been active since ∼ 61 ka, the difference in the geological slip rates derived from the difference of the total slip (∼ 600 m) could be as high as ∼ 10 mm yr⁻¹, which is more than half of the total shortening rate across the central Himalaya measured geodetically (Fig. 6). This differential rate could be even larger if the structure is younger or less if the structure is older. However, even if the fault had been activated 5 times earlier (∼ 300 ka), the differential rate would still account for ∼ 10 % of the overall geodetic shortening (Fig. 6). The differences in estimates of the total slip thus have a knock-on effect on seismic hazard assessments: the overestimation of the total slip along the Bardibas thrust leads to underestimation of the slip rate deficit of the MFT and thus underestimated seismic potential. Moreover, scenarios of slip partitioning among faults and the implications for structural evolutions may vary. For example, the discrepancy between the slip of the Patu thrust, a hinterland strand of the MFT, and the Bardibas thrust may actually be more significant than 1.8 km for the former and 1.7 km for the latter estimated by Almeida et al. (2018).

Except the differential estimates of the total slip, another source of uncertainty for estimating geomorphic slip rates by terraces above is the solely decided ramp dip. Almeida et al. (2018)'s seismic image in their Fig. 2 is described as providing direct evidence of fault dip, which is claimed to allow for a more precise conversion between uplift and input displacement. Yet, the seismic image allows alternative interpretations, and if the ramp dip is 10^∘ steeper or gentler than the drawn dip of ∼ 20^∘, then the ratio of crestal uplift to the input slip could differ by ∼ 30 %–60 % for all types of fault-bend folds (Hardy and Connors, 2006). However, the potential uncertainties of the fault geometry inherited from the seismic image are neither expressed nor included in their interpretations, so future applications of Almeida et al. (2018)'s results may inherit (unawares) these uncertainties.

4 Discussion and conclusion

In this study, I resolve the total geological slip of a shear fault-bend fold with a simple solution from Suppe et al. (2004), validated with a case from Le Béon et al. (2019). The benefit of the proposed solution is that the full geometry does not need to be visualized; instead, the solution relies on the ramp dip and features (i.e. backlimb length and dip), which are commonly located in the better-imaged shallow part of seismic images. Though this solution may entail further considerations in more structurally complicated cases (e.g. structures with changing bed lengths, a significant amount of off-fault deformation, or curved fault planes) and the degrees of uncertainty may vary from case to case (e.g. Randle et al., 2018), I quantitatively demonstrate a previously neglected source of uncertainty in estimating geological slip rates with a simple case of the Bardibas thrust in the central Himalaya. This subtle difference in structural interpretations and estimates of the total slip (∼ 600 m) for the structure can substantially affect the following estimates of the geological slip rate by up to ∼ 10 mm yr⁻¹ and thus implications for seismic potential and structural evolutions, especially if the structure is younger. Also, the exclusively interpreted ramp geometry could lead to up to ∼ 60 % difference in the uplift-to-slip ratio. Therefore, I suggest that seismic reflection can provide more precise structural geometry and estimates of the total slip only if the method for interpretations is properly applied and compatible with available data quality. Otherwise, it may introduce unintended uncertainty.

Structural geometry is central to estimating geological or geomorphic slip rates and thus to assessment of seismic hazard. The unclear fault geometry at depth (Hu and Stevens, 2022) or near the surface, as shown in this study, can lead to non-negligible uncertainty in terms of the issues mentioned above. Though ambiguous fault geometry is unavoidable, subjectivity in structural interpretations can be reduced by either including uncertainties in interpretations or proposing multiple working hypotheses (Chamberlin, 1965).