Inverse trishear modeling of bedding dip data using Markov chain Monte Carlo methods

Journal of Structural Geology 80 (2015) 157e172

Contents lists available at ScienceDirect

Journal of Structural Geology journal homepage: www.elsevier.com/locate/jsg

Inverse trishear modeling of bedding dip data using Markov chain Monte Carlo methods David O.S. Oakley*, Donald M. Fisher Department of Geosciences, Pennsylvania State University, University Park, 16802, PA, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 30 April 2015 Received in revised form 14 September 2015 Accepted 21 September 2015 Available online 26 September 2015

We present a method for fitting trishear models to surface profile data, by restoring bedding dip data and inverting for model parameters using a Markov chain Monte Carlo method. Trishear is a widely-used kinematic model for fault-propagation folds. It lacks an analytic solution, but a variety of data inversion techniques can be used to fit trishear models to data. Where the geometry of an entire folded bed is known, models can be tested by restoring the bed to its pre-folding orientation. When data include bedding attitudes, however, previous approaches have relied on computationally-intensive forward modeling. This paper presents an equation for the rate of change of dip in the trishear zone, which can be used to restore dips directly to their pre-folding values. The resulting error can be used to calculate a probability for each model, which allows solution by Markov chain Monte Carlo methods and inversion of datasets that combine dips and contact locations. These methods are tested using synthetic and real datasets. Results are used to approximate multimodal probability density functions and to estimate uncertainty in model parameters. The relative value of dips and contacts in constraining parameters and the effects of uncertainty in the data are investigated. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Trishear Fault-propagation folds Fold kinematics Markov chain Monte Carlo

1. Introduction Trishear is a kinematic model for fault-propagation folds, in which a triangular zone of distributed deformation occurs ahead of the fault tip. It is capable of reproducing features such as nonuniform forelimb dips and footwall synclines (Erslev, 1991) d features that are not explained by kink band models (Suppe, 1985; Suppe and Medwedeff, 1990) but are frequently observed in the field. Deformation within the trishear zone can be described in terms of a velocity field (Hardy and Ford, 1997; Zehnder and Allmendinger, 2000) in which the velocity of a particle at any given point can be calculated from its position relative to the fault tip. By integrating this velocity through the slip of the fault as the fault tip propagates, the geometry of a trishear fold can be modeled from the initial stratigraphy and the trishear model parameters: tip position, fault dip, fault slip, trishear zone apical angle, propagation to slip ratio, and the concentration factor, which if greater than one concentrates deformation toward the center of the trishear zone. No analytic solution is known, but the integration can be performed

* Corresponding author. E-mail addresses: [email protected] (D.O.S. Oakley), [email protected] (D.M. Fisher). http://dx.doi.org/10.1016/j.jsg.2015.09.005 0191-8141/© 2015 Elsevier Ltd. All rights reserved.

numerically. Balanced cross sections, although fundamental to structural geology, are often non-unique, and multiple interpretations may be possible. This is especially likely when subsurface data are absent. Interpretations based on trishear kinematics are by no means exempt from these concerns, and the six or seven (if the concentration factor is included) free parameters in the model make it likely that there will be a range of reasonable interpretations. Interpretations, such as shortening estimates based on a balanced cross section, may have substantial uncertainty (Judge and Allmendinger, 2011) and will be more useful if this uncertainty can be quantified. Trishear folds, which must be modeled numerically, are amenable to methods that require testing a large number of potential cross-sections. Such techniques have been used to find a best-fit trishear model for given data beginning with the work of Allmendinger (1998), who used a grid search over the parameter space. The grid search approach has subsequently been employed by other authors (e.g. Allmendinger and Shaw, 2000; Allmendinger et al., 2004; Cardozo, 2005; Lin et al., 2007), principally with the aim of identifying a best-fit model. Quantifying uncertainty in trishear parameters presents a greater challenge than finding a best-fit model alone, but it is essential to a full understanding of a structure. Although not

158

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

commonly done, formal uncertainty can be estimated from grid search results, as we demonstrate below. Other methods to quantify uncertainty in trishear model parameters have been proposed by Cardozo and Aanonsen (2009) and Regalla et al. (2010). Cardozo and Aanonsen (2009) use the randomized maximum likelihood (RML) method, in which the best fit is found by optimization for many different realizations of the data and the results are plotted as histograms for each trishear parameter. Regalla et al. (2010) use a Monte Carlo simulation and plot histograms of all models for which the objective function is below a certain threshold. Cardozo et al. (2011) have also shown that simulated annealing can be used to estimate the range of possible models. In this study, we propose the use of Markov chain Monte Carlo methods. Such methods, while not previously applied to trishear, are commonly used for data inversion and are capable of both finding a best fit and estimating uncertainty (Tarantola, 2005). Trishear kinematics is fully reversible and so can be used to test possible cross-sections by restoration (Allmendinger, 1998). A bed trace, imaged in seismic data or exposed in outcrop along the length of the fold, can be restored in this manner, and the model results can be compared to a predicted geometry using an objective function (such as distance from a point to a best fit straight line) in order to evaluate the goodness of fit. Mapped contacts along a surface profile can also be restored and matched across the fold or to a known original depth, but dip measurements pose more of a problem. The approach used in previous work to match dip measurements to trishear models (Cardozo, 2005; Regalla et al., 2010) has been to forward model beds, interpolate dips between modeled points, and attempt to match the observed dips. In summary, work by various authors has demonstrated the value of the inverse method, using a variety of inversion algorithms, in fitting trishear models to data. More limited, but significant, work has shown that the technique can be applied to datasets consisting of dips or of mixed dip and point data and has proposed some methods by which error in trishear parameters can be estimated. We build on this body of knowledge by developing a method for the direct restoration of dip data, testing Markov Chain Monte Carlo techniques for trishear data inversion and error estimation, and investigating the relative value of contact positions and bedding dips in constraining a trishear model.


   v0 jyj 1 s þ 1 sgnðyÞ mx 2  
 v0 m jyj ð1þsÞ s vy ¼ 1 2ð1 þ sÞ mx vx ¼

xm  y  xm; s  1 where sgn(y) indicates the sign of y, and s is a concentration factor. Increasing s concentrates deformation toward the center of the trishear zone. When s ¼ 1, the field is termed “linear” (Zehnder and Allmendinger, 2000) or “homogeneous trishear” (Erslev, 1991). This velocity field properly describes the velocities in the (z,h) coordinate system. In the (x,y) coordinate system there is an additional component of relative motion in the x direction, equal to ev0(P/S), due to fault propagation. P/S is the ratio of fault propagation to fault slip. In summary, there are seven parameters that determine the form of a trishear fold: the two coordinates of the fault tip (xt, yt), the fault ramp angle (fault dip), the total slip on the fault, P/S, f, and s. This velocity field can be used directly to restore a folded bed by incrementally moving points along the bed. For dip data, the forward modeling approach (Cardozo, 2005; Regalla et al., 2010), requires moving a large number of points for each dip, which is timeconsuming. If dips are instead restored directly, via an equation for the rate of change of dip, dip data can be inverted in the same manner as point data. This approach is much more efficient, allows a restoration approach to be used with datasets that contain both dip and point (bed or contact) data, and also reduces the errors inherent in interpolating between points. The velocity field of Eq. (1), or any trishear velocity field, can be used to define a strain rate tensor in two-dimensions:

3 vvx vy 7 7 7 vvy 5

2

vvx 6 vx 6 e_ ¼ 6 4 vvy vx

(2)

vy

If a dip in cross-section is defined by an angle that it makes with the x-axis in the trishear coordinate system, q, then the strain rate tensor can be transformed into a coordinate system aligned with the dip: 0

_ e_ ¼ aT ea 2. Methods 2.1. Velocity equations The trishear velocity field proposed by Zehnder and Allmendinger (2000), while not the only possible velocity field, is among the simplest and most widely used, often in its linear form (s ¼ 1 in Eq. (1)) (Hardy and Allmendinger, 2011). This formulation is supported by the finite element modeling of Cardozo et al. (2003) for an elastoplastic material, which produces a similar fold shape. This velocity field is defined with reference to a coordinate system (x,y), for which the origin is at the fault tip and moves with it as the fault propagates. A second coordinate system (z,h) has its origin fixed at the initial fault tip position. The trishear zone is a triangular region defined by the angle (4) between its boundaries and the xaxis. The tangent of 4 is denoted m. In some cases, the trishear zone may be asymmetric, requiring two 4 values (Zehnder and Allmendinger, 2000), but we will focus on the symmetric case. The hanging wall velocity outside the trishear zone is labeled v0. Fig. 1 shows the geometry of the trishear zone with the coordinate axes and important variables indicated. The trishear velocity field of Zehnder and Allmendinger (2000) is:

(1)

(3)

_ is the strain rate tensor in a Cartesian coordinate system where e′ with its x0 -axis parallel to the dip, and a is the rotation matrix:

 a¼

cosq sin q

sinq cosq

 (4)

_ 21 of the rotated strain rate tensor will The off diagonal term e′ then be the rate of rotation of a line parallel to the x0 -axis (Allmendinger et al., 2012) and will therefore be equal to the rate of change of dip. 0 vvy vvy vvx vvx 2 cosqsinq þ cos2 q  cosqsinq q_ ¼ e_ 21 ¼  sin q þ vx vx vvy vy

(5) For the trishear velocity field of Eq. (1), this simplifies to

v q_ ¼  0 2sx

" #2 rffiffiffiffiffi
1þs 1s
pffiffiffiffiffi jyj 2s 1 jyj 2s m cosq  sgnðyÞ sinq m mx mx

(6)

Note that Eq. (6) uses the convention that q is positive

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

159

Fig. 1. The geometry of the trishear zone, trishear coordinate system, and model parameters.

counterclockwise up from the x-axis, meaning that in cross-section with x positive to the right, dips down to the left will be positive, while those down to the right will be negative. This convention can be readily switched if desired. Eq. (6) can be used along with Eq. (1) to incrementally move and rotate a point, thus allowing a dip measurement to be restored to its pre-deformational position and dip. 2.2. Uncertainty analysis For any inversion of bed or dip data, each model test produces a measure of misfit according to the objective function. This misfit can be translated into a probability, P(djp), that is the probability of obtaining the observed data (d) given the parameters of the tested model (p). If errors in the data are uncorrelated and normally distributed with uncertainty s, this probability will be:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u
1 !2
1þs !2 u v0 Dt mv0 Dt jyj s jyj s t 2 sgnðyÞ dytþDt ¼ 1þ dyt þ dx2t 2sx 2sx mx mx Since jyj  mx within the trishear zone, the change in dx and dy with each time step will usually be small, especially for large values of x. The uncertainty in the dip, dq at time step t þ Dt will be a function of the previous dq and also of uncertainties in the position (dx, dy), since q_ is a function of q, x, and y. The equation for error propagation for q is therefore

dqtþDt

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u vq_ vq_ vq_ t ¼ 1 þ Dt dqt þ Dt dxt þ Dt dyt vq vx vy

(10)

where c2

PðdjpÞfe 2

(7)

where c2 is the sum of squares of the errors divided by s2. Assuming no a priori information and ignoring uncertainties in the theoretical relationship between d and p, P(djp) will be proportional to P(pjd) or the probability of the model parameters given the data (Tarantola and Valette, 1982; their Eq. (6.9)). If one does have prior information, this can be multiplied in to the final probability. Uncertainties are typically estimated for data in the deformed state, but the inverse (cross-section restoration) method produces an error for the undeformed state. It is therefore necessary to propagate errors in position or dip through deformation in the trishear zone. If trishear deformation is solved numerically in a series of timesteps of length Dt or slip increments of length Dx ¼ v0*Dt, then the uncertainties after each step can be derived from Eq. (1) for points moving through the trishear zone and from Eq. (6) for dips. For points, the uncertainties in x and y, here termed dx and dy, are:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u
1 !2
1s !2 u v Dt v0 Dt jyj s jyj s t dxtþDt ¼ 1  0 sgnðyÞ dx2t þ dy2t mx 2sx 2msx mx and

(8)

 
1 
vq_ v jyj s 1 y x ¼ 0 sgnðyÞ  sinð2qÞ þ cosð2qÞ ; vq sx 2 x y mx
1=s  y vq_ v jyj ¼ 20 2 sgnðyÞ cos2 q ð1 þ 2sÞ vx 2s x x mx  x  ð2 þ 2sÞsin q cos q þ sin2 q y

(11)

(12)

and


1 "
 vq_ v0 x jyj s ¼  2 2 sgnðyÞ sin q cos q ð1 þ sÞcos2 q  2 vy y mx 2s x #
2 x þ ð1  sÞ sin2 q y (13) where, as in Eq. (6), q is considered positive for dips down to the left in the cross section and negative for dips down to the right. If the uncertainties in position of dip measurements are considered to be negligible relative to the uncertainty in dip, then

160

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

including grid search and multiple forms of Markov Chain Monte Carlo simulation. We apply this program first to a synthetic section and then to a natural example.

Eq. (10) simplifies to

dqtþDt ¼ 1 þ


1 
 v0 Dt jyj s 1 y x sgnðyÞ  sinð2qÞ mx sx 2 x y !

þ cosð2qÞ

(14)

dqt

3. Results 3.1. Synthetic section

This might occur if the positions of dip measurements are well constrained by GPS but more significant uncertainty exists in the dips themselves. It is also suggested by the fact that Eq. (11) has an x1 term where Eqs. (12) and (13) have x2, which will make the effect of uncertainty in position much smaller than that of uncertainty in dip for large values of x. For this paper, we use this simplification. This interpretation of the error in the restored model in terms of probability means that each model tested is sampled from a density function over the parameter space, which when normalized is the probability density function for the parameters. Given that this density function exists in six dimensions when solving for fault tip position, total slip, ramp angle, P/S, and f, it cannot be visualized directly, but by integrating over this density function, one can calculate the marginal probability density functions for each parameter as well as statistics such as expected value, standard deviation, and covariance (Tarantola and Valette, 1982). The probability density function allows the uncertainty and range of possible values for each parameter to be characterized, whether or not they follow a Gaussian distribution. With a grid search algorithm, one can readily integrate over the grid of probabilities using quadrature. Alternatively, the treatment of the error as a probability lends itself to solution by Markov Chain Monte Carlo methods, which sample directly from the final probability distribution. To test these new techniques, we have developed a program for trishear inverse modeling, using Eqs. (1) and (6) to restore points and dips respectively and Eqs. (8), (9) and (14) to propagate errors. The program is capable of using several inversion techniques,

To test the techniques presented here, we first use data derived from a synthetic cross-section. This approach allows us to test the algorithms in a controlled case before applying them to natural, and thus more uncertain, structures. The synthetic section (Fig. 2) was created using the FaultFoldForward v. 6.1.1 program by Allmendinger (http://www.geo.cornell.edu/geology/faculty/RWA/ programs/faultfoldforward.html), so that the program used for inversion could be validated against this older, more-tested forward-modeling trishear program. Length is in generic units, which could be scaled as appropriate to be representative of folds from outcrop to map scale. Beds were initially horizontal at elevations from 50 to 250 units, spaced every 10 units. Modeled points were spaced every 1 unit along the beds, and fault slip proceeded in 1 unit increments. Model parameters are shown in Fig. 2. Data were extracted from the model along a line at a 10 slope, as shown in Fig. 2. This line was chosen as one that would intersect all 21 modeled beds and pass from the hanging wall into the trishear zone. Contact locations, where the line intersected a bed, were taken from every fifth bed starting with the lowest, and dips were calculated at every bed, using the slope between the two nearest modeled points. In this manner, a simulated profile of dip and contact data was produced. The simulated profile was then perturbed to simulate measurement error, using a Gaussian distribution with standard deviation of 5 for dips and 5 units (horizontally and vertically) for contact positions. The positions of dip measurements were not altered. Inversion was first performed using a grid search. Grid

Fig. 2. The synthetic cross section and the model parameters used to create it. The bold line is the profile along which data were taken.

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172 Table 1 Grid used for the grid search algorithm. Parameter

Minimum

Maximum

Step size

xt yt Total Slip Ramp Angle

300 50 0 10 10 0

700 250 400 50 50 3

10 10 10 2 2 0.1

f P/S

parameters are shown in Table 1. The concentration factor, s, was assumed to be 1, and a model was generated to fit for all six other parameters. In total, 482,599,971 models were tested. This took 56,247 s (about 15.6 h), running on a desktop computer with an Intel Core i7-4790K four-core, 4 GHz processor. All times quoted in this paper are for this same computer. For the grid search, eight processes were run in parallel using the OpenMP API for parallel processing, for a total processor time of about 125 h. The undeformed stratigraphy was assumed to be flat with known bed elevations, so dip errors were simply the restored dips, and errors for the contacts were the distance between each restored contact point and a line representing the expected position of the corresponding horizon, which for initially horizontal stratigraphy is the difference in elevation between the point and the horizon. Probabilities were calculated using uncertainties in the data of 5 for dips and 5 units for contacts, which were propagated as described above. The results were integrated to calculate probability densities (Fig. 3) and statistics (Table 2). Several important results should be noted. First, despite the large number of models run, several of the marginal probability distributions come to sharp peaks, suggesting that the grid is still not fine enough to accurately capture the curvature of the distribution. Second, several of the distributions are multimodal, having peaks not only at the correct solutions, but also at alternative solutions. Since means and standard deviations of the multimodal distribution are not very useful in estimating model parameters, we also list these statistics for the two groups of solutions, dividing

161

them into those for which P/S  1 and those for which P/S  1. Given the difficulties with the grid search approach, there is a need for a faster approach that is still capable of accurately estimating the probability density functions for the model parameters. Markov chain Monte Carlo methods are frequently used in data inversion to achieve this purpose and are suited to the probabilistic treatment of the model errors described above. One difficulty with the basic MetropoliseHastings algorithm (Metropolis et al., 1953; Hastings, 1970) in this case is that without knowing the uncertainty in the final parameter values, it can be difficult to choose a good proposal distribution from which the random walk of the algorithm will choose new model parameters. To overcome this difficulty, we first used the Robust Adaptive Metropolis (RAM) algorithm (Vihola, 2012), which adapts the proposal distribution in order to achieve an ideal acceptance rate for proposed models in the Markov chain. Histograms of the results of the RAM algorithm are shown in Fig. 4 for a run consisting of 1 million models, which took only 665 s (or ~11 min). The initial model was randomly chosen from within the parameter space, which has the same maximum and minimum values as the grid search parameter space. Dashed lines in Fig. 4 show the correct values of model parameters that were used to make the original model. Model peaks lie near these correct values. Where the RAM algorithm falls short is that the second set of peaks, as seen in Fig. 3, are not observed here. Thus the algorithm has not captured the full range of models that can fit the data. To reproduce the multimodal distributions seen in Fig. 3, we made use of the adaptive parallel tempering (APT) algorithm of Miasojedow et al. (2013). This method is another type of Markov chain Monte Carlo algorithm, in which multiple chains are run targeting increasingly higher “temperature” versions of the probability distribution, meaning that the higher level chains move more easily through the parameter space. Chains are allowed to swap states, thus allowing the lower level chains to explore multiple probability maxima, and both the temperature distribution and the proposal distribution for each level are adapted to achieve an ideal acceptance rate. For this, we follow Vihola (2012) and

Fig. 3. Probability density functions for synthetic model parameters generated from a grid search. The dashed lines show the parameter values used to create the original model.

162

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

Table 2 Means (m) and standard deviations (s) calculated from the grid search results for the synthetic section. Parameter

xt yt Total Slip Ramp Angle

f P/S

P/S  1

Full Model

P/S  1

m

s

m

s

m

s

531 131 195 32.4 32.7 1.05

39.6 25.6 20.9 3.8 2.9 0.54

483 103 196 32.2 33.0 0.41

9.0 8.1 20.3 3.8 3.3 0.13

562 150 194 32.5 32.6 1.48

10.4 12.7 21.3 3.8 2.6 0.099

Miasojedow et al. (2013) in using a target acceptance rate of 23.4% for both parts of the algorithm. Adaptation of the proposal distributions is performed using the RAM method. The results of the lowest level, which is untempered, are taken to represent the probability density of the model parameters. Histograms of the results of the APT algorithm are shown in Fig. 5a and statistics of the distributions are given in Table 3 for a run consisting of 8 energy levels, each with 1 million models, starting with a randomly chosen initial model. To run these models took 1224 s (~20 min), or longer than the RAM method but much less time than the grid search. Parallel processing was used to run models at different energy levels at the same time, but unlike the grid search, this algorithm requires frequent returns to sequential processing. Fig. 5a and Table 3 show that this method is able to reproduce the multimodal distributions seen in the grid search results (Fig. 3 and Table 2), but it does so with far fewer models and with a higher resolution in areas of interest. The lowest energy chain switches regularly between the two modes (Fig. 5b), allowing both to be accurately represented in the histograms. A comparison of the original model used to produce the synthetic section with the best-fitting P/S < 1 model (Fig. 6) shows that both solutions are similar. The fit is best near the surface profile line, where it is constrained by data. If a different line is used, the same effect is

seen: alternative models match the original model well along the chosen line but may differ from it away from that line. The difference between the two models in Fig. 6 is greatest in the forelimbs of the folds, where the curvature is greatest and where the shape of the fold is most sensitive to the trishear parameters. Good data in the forelimb are thus necessary to distinguish between possible models. Since trishear parameters may covary, contours of twodimensional histograms are plotted in Fig. 7. The division between models with P/S < 1 and P/S > 1 appears clearly in plots including P/S and in the yt vs. xt plot, where the two sets of models define two distinct regions in which the fault tip may be found. There is a positive correlation between xt and yt, showing linear regions in the cross-section space where the tip is most likely to be found. There is also a strong negative correlation between ramp angle and total slip, which is likely a product of the need to restore the contacts by a set vertical distance to their undeformed elevations. An approximately inverse correlation between P/S and f has been noted by previous authors (Allmendinger et al., 2004; Cardozo, 2005; Hardy and Allmendinger, 2011). In Fig. 7, we observe this behavior when P/S > 1, but we see a positive correlation between the two variables when P/S < 1. Thus, more generally, P/S approaches 1 as f increases. The use of both dips and contacts along the surface profile provides more data than either alone. APT model runs of 1 million models each were produced using each type of data separately. When using just the points for the five contacts, the probability distributions are much wider, indicating a loss of precision (Fig. 8a). Peaks remain near the correct values, but may be offset somewhat, as is particularly noticeable for total slip. The bimodal nature of the distribution remains visible in the xt and P/S histograms, but the two peaks are less distinct from each other. When using just dip data, precision is also lost and the losses in accuracy are more severe (Fig. 8b). The histogram for f, for instance, has its mode at about 45 and remains high at the boundary of the parameter space, but it is low at the correct value of 30 .

Fig. 4. Histograms of results from the RAM algorithm. Each histogram is divided into 100 equally spaced bins. The total number of models is 1 million. The dashed lines show the parameter values used to create the original model.

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

163

Fig. 5. (a) Histograms of results from the APT algorithm. Each histogram is divided into 100 equally spaced bins. The total number of models is 1 million. The dashed lines show the parameter values used to create the original model. (b) Plot of the course taken by the lowest energy chain of the APT algorithm.

To compare the effect of propagating the uncertainty to that of assigning an uncertainty to the restored data, we ran the APT algorithm for the bed and dip data individually without propagating uncertainty. Uncertainties of 5 and 5 length units were then used to calculate a probability from the errors in the restored data. Results are shown in Fig. 9a and Fig. 9b. For contacts, the histograms in Fig. 9a are not very different from those in Fig. 8a, indicating that the effects of propagating errors are minimal. As noted above, Eqs. (8) and (9) suggest that the change in uncertainty is likely to be

small with each slip increment, and these results support that assumption. For dips, however, this assumption is not valid. Histograms for total slip, ramp angle, and f in Fig. 9b all have peaks at the limit of the allowed parameter space for high values of slip and low values of ramp angle and f. These results are not close to the correct values that were used to make the model, and are considerably less accurate than the results obtained with error propagation (Fig. 8b). The data uncertainty will have an important effect on the final

164

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

Table 3 Means (m) and standard deviations (s) calculated from the APT results for the synthetic section. Parameter

xt yt Total Slip Ramp Angle

f P/S

P/S < 1

Full Model

P/S > 1

m

s

m

s

m

s

529 130 193 32.7 32.3 1.05

38.1 24.2 20.8 4.0 3.2 0.54

484 103 193 32.8 32.6 0.400

8.9 7.8 20.0 3.9 3.1 0.13

559 148 194 32.6 32.1 1.48

8.7 10.3 21.3 4.1 3.2 0.10

probability distribution. In this example, the uncertainties of 5 for dips and 5 units for contact positions were those used to perturb the data in the first place, but in a natural example they will need to be estimated based on the quality of the measurements. Fig. 10a and Fig. 10b show the results, from the APT algorithm, when uncertainties of half and twice these values respectively are used instead. As expected, the uncertainty in model parameters increases and decreases with the uncertainty in the data. When using half the original uncertainties, the difference between the two peaks in the histograms is greater than in the original model, with the correct value at P/S z 1.5 clearly preferred over the P/S z 0.5 solution. When the uncertainties are twice their original values, on the other hand, the two peaks, in addition to being broader (less precise), are more nearly equal in amplitude. Thus differences in the uncertainty of the data will affect not only the precision of the results, but also the degree to which alternative models appear likely. 3.2. Waterpocket monocline Cardozo (2005) uses the Waterpocket Monocline, a Laramide structure in southern Utah, as an example in the paper that first introduced a method for fitting trishear models to surface profile data, including dip data. As we are proposing new methods for this

purpose, we apply our methods to this same example for comparison. The data are taken from Bump (2003) and consist of dips and formation contacts mapped along a topographic profile as well as contact points on a footwall stratigraphic column, derived from well logs and projection from the surface. We assume that the undeformed stratigraphy is flat, with the bed contacts at the elevations shown in the footwall stratigraphic column of Bump (2003). While Cardozo's (2005) method used only the dips as measured at the contacts, we use all the dips shown in the profile, since by our method they can be restored individually. We also use a larger parameter space (Table 4) than that used in Cardozo's (2005) grid search. We use the APT algorithm with 100,000 models and 10 energy levels. This run took about 57 min; this was significantly slower than the synthetic model, due to the much larger slip distances, and thus larger number of steps for deformation within the trishear zone. The first 20,000 models were removed as a burn-in period before plotting histograms, since this was approximately the number needed for the Markov chain to find the regions of high probability after starting from a randomly chosen value within the large parameter space. Uncertainties are 5 for dips and 25 m for contacts. 5 is an estimate of the likely uncertainty in dip measurements, while 25 m is an approximation of the uncertainty in our digitization of contact positions from the surface profile in Cardozo (2005). Fig. 11a shows histograms of the inversion results for the Waterpocket monocline. Like those for the synthetic section, these results show two significant modes, one with P/S > 1 and one with P/S < 1. Unlike the synthetic section inversion, in which the Markov chain switched frequently between the two solutions (Fig. 5b), in this case the chain switched from the P/S > 1 solution to the P/S < 1 solution only once, due to the great distance between the two (Fig. 11b). A further test (with one million models) showed that after this switch, the chain did not switch back. Thus the relative magnitude of the two peaks is not in this case a good indicator of their relative probability; it is only the fact that both exist and provide good solutions that should be noted. In fact, the best

Fig. 6. Comparison of original synthetic model (solid lines) and best-fitting P/S < 1 solution (dashed lines).

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

165

Fig. 7. Contours of two-dimensional histograms of APT results. The lowest contour in each plot is 1000 models. Contour intervals, clockwise from top left, are: 4000, 2000, 2000, and 1000 models. The full parameter space was divided into 100 bins in each dimension, but only the regions of maximum probability are shown here.

models in the P/S < 1 group produced somewhat higher probabilities than the best P/S > 1 models. The P/S < 1 solutions correspond to ramp angles very near vertical. Models in this category have an average ramp angle of 89.9 . While mathematically valid as a solution for trishear kinematics, these solutions are mechanically unrealistic. Reverse slip on an essentially vertical fault is inconsistent with the subhorizontal direction of maximum principal stress inferred for the Laramide orogeny (Erslev and Koenig, 2009). This geometry is also inconsistent with other Laramide faults, which seismic and well data show to dip less steeply (e.g. Smithson et al., 1979; Stone, 1993), consistent with horizontal shortening. Results with P/S > 1 are relatively similar to the best fit model of Cardozo (2005). The means and 1s uncertainties for the parameters (xt, yt, slip, ramp angle, f, P/S) are: (6637 ± 11.9 m, 1682 ± 19.9 m, 3157 ± 78.9 m, 43.9 ± 1.4 , 58.0 ± 1.3 , 2.36 ± 0.013), and the bestfitting model is (6638 m, 1688 m, 3173 m, 43.7, 57.8 , 2.37). For comparison, Cardozo's (2005) grid search best fit is (6600 m, 1270 m, 3800 m, 35 , 52.5 , 2.25). These values are outside the uncertainty range estimated for our parameters but still within the same general area of the parameter space. Given the different inversion methods (APT vs. grid search), different objective functions (probability vs. chi-square statistic), differences in the parameter space (our yt and ramp angle means are outside the range tested by Cardozo (2005)) and use of additional dip data from the Bump (2003) profile, it is not surprising that there is some difference between the results. Fig. 12 shows a comparison of the folds produced by our best-fit model and that of Cardozo (2005), plotted on top of the surface profile from Cardozo (2005).

In the Waterpocket example, and in many other natural cases, assigning an uncertainty to the data may be difficult. In addition to measurement errors, which may or may not be well known, the complexities of a natural structure mean that it is unlikely to perfectly follow the trishear model. To explore the effects of different uncertainty values on the results, we ran APT inversions for uncertainties in dip of 2.5 , 5 , and 10 and uncertainties in contact position of 12.5 m, 25 m, and 50 m, for a total of 9 different combinations. Each inversion used 100,000 models and 8 energy levels. Ramp angle was limited to values less than 70 , in order to have only a single peak in which to measure uncertainty, and the rest of the parameter space was as in Table 4. The first 50,000 models were removed as burn-in time, as some runs took nearly this long to reach the region of interest, leaving the last 50,000 models from each run to be used to calculate the means and standard deviation. Means are shown in Table 5 and standard deviations in Table 6. As expected, standard deviation increases as uncertainty increases, but the means are not very sensitive to the change. In many cases, the means for the other uncertainty values are within the 1 s uncertainty of the means for the 5 and 25 m uncertainties, and in no cases do they jump to a wholly new part of the parameter space. 4. Discussion The restoration approach for dip data and the APT inversion algorithm show substantial promise in solving trishear problems. The superiority of APT over the grid search is straightforward. It can find the probability distributions of model parameters with higher

166

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

Fig. 8. Histograms of APT results fit to (a) surface contacts only and (b) dip data only. The total number of models is 1 million. The dashed lines show the parameter values used to create the original model.

resolution than the grid search and with many fewer models and thus less time required. The large number of models used in the grid search for Fig. 3 would not be possible if the number of data or slip distance were significantly larger, since these are the two main factors controlling how many times Eqs. (1) and (6) must be evaluated. Even with this many models, the grid search produces far less resolution in the areas of interest in the probability density functions than does the APT algorithm. If one were searching for

additional parameters beyond the six in this example, the grid search would become even more impractical. Such additional parameters might include s, depth to detachment, or geometry of a listric fault as in Cardozo and Brandenburg (2014). APT can identify both local and global probability maxima, but in comparison to optimization methods, it has the advantage of also being able to determine uncertainty in model parameters, and it is less likely to be susceptible to the local minima that those methods can be

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

167

Fig. 9. Histograms of APT results fit to (a) surface contacts only and (b) dip data only, without propagating errors. The total number of models is 1 million. The dashed lines show the parameter values used to create the original model.

trapped by Cardozo and Aanonsen (2009). Nonetheless, if only a best fit is desired and the problem is well constrained, an optimization method will likely be faster. Other methods may also be chosen depending on the situation. For example, Cardozo and Aanonsen (2009) successfully used a randomized maximum likelihood method to estimate uncertainties in four parameters for a trishear model of the Santa Fe Springs anticline based on seismic reflection data. Markov chain Monte Carlo techniques that are

simpler, and faster, than APT may also be useful in many situations. The results shown in Fig. 4 suggest that the RAM method would be sufficient if the parameter space were sufficiently constrained that a multimodal distribution was not likely. Like these other approaches, the APT algorithm can also be used to find best-fit models and to estimate uncertainty in well-constrained problems, but it is particularly suited to exploring a large, poorly constrained parameter space and to reproducing multimodal probability distributions.

168

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

Fig. 10. Histograms of APT results. (a) s ¼ 2.5 for dips and 2.5 length units for contacts. (b) s ¼ 10 for dips and 10 length units for contacts. The total number of models is 1 million. The dashed lines show the parameter values used to create the original model.

The multimodal nature of the probability density function for our synthetic model is important to note. If it were not known a priori which solution is correct, as would be the case for a natural example, it would be necessary to consider both possible solutions as valid. Such results might indicate the need for additional data in order to distinguish between the two models. In this case, the two models are divided most clearly by their P/S values. In published studies of real structures, most P/S values are greater than 1 (Pei

et al., 2014), and it may be that values less than 1, while mathematically possible are uncommon in nature. More studies of natural trishear folds and a better understanding of the physical underpinnings of the propagation-to-slip ratio are needed to show whether or not this is the case. At present, however, care must be taken when assuming a particular trishear model for a fold if more than one model is possible. Bedding and contact locations along a surface profile work best

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172 Table 4 Parameter space for the Waterpocket Monocline inversion. Parameter

Minimum

Maximum

xt yt Total Slip Ramp Angle

5000 1000 2000 0 0 0

8000 2000 8000 90 90 5

f P/S

169

when used together to constrain a trishear model. Dips alone produce a less accurate result (Fig. 8b) and either data type alone results in significantly lower precision. Different data can help to constrain different parameters, as for example dips alone provide a much more precise constraint on P/S (Fig. 8b) than do contacts alone (Fig. 8a), while as noted above, dips provide a poor constraint on f. The inaccuracy in fitting dip data alone is much greater without error propagation (Fig. 9b) than with (Fig. 8b). The results shown in

Fig. 11. (a) Histograms of APT model results for the Waterpocket Monocline. The total number of models is 80,000. Histograms are divided into 100 bins. (b) Plot of the course taken by the lowest energy chain of the APT algorithm.

170

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

Fig. 12. Modeled cross sections of the Waterpocket Monocline, using our best fit model (solid line) and Cardozo's (2005) best fit model (dashed line). The background is from Fig. 4 of Cardozo (2005), which is in turn derived from Bump (2003).

Fig. 9b are a more direct representation of the error in the restored dips than are those shown in Fig. 8b. This is because in Fig. 9b the dip errors are converted into probabilities with the same s for all models. Thus this approach is analogous to methods commonly used for restoring beds (Allmendinger, 1998) that seek to minimize the error in the restored beds. The maximum probability model in Fig. 9b is also the model with the lowest RMS error in the restored dips. Specifically, RMS error for the restored dips for this model is 1.6 , while the model parameters (xt, yt, slip, ramp angle, f, P/S) are (501, 187, 392, 11.1, 10.5, 0.79). In contrast, the best-fit model in Fig. 8b has a restored dips RMS error of 9.9 with parameters of (574, 173, 231, 46.5, 45.0, 1.28). Although not perfect, these

parameters are closer to the correct values shown in Fig. 2, despite the much higher RMS error in the restored dips. For comparison, the correct parameter values give an RMS error of 4.1 in restored dips, due to the fact that the original dips were perturbed to simulate measurement errors. For beds or contacts, the lack of major differences between Figs. 8a and 9a indicates that error propagation is not critical and minimizing RMS error in the restored beds is appropriate. For dips, however, error propagation is necessary. While both show significant losses of accuracy and precision over an inversion with both dip and contact data, the results with error propagation (Fig. 8b) are clearly more accurate than those without (Fig. 9b). Therefore, simply minimizing RMS

Table 5 Mean values of model parameters for different uncertainties in dips (sdip) and contact location points (spoint).

Table 6 Standard deviations for model parameters for different uncertainties in dips (sdip) and contact location points (spoint).

spoint

sdip 2.5

12.5 25 50 12.5 25 50 12.5 25 50 12.5 25 50 12.5 25 50 12.5 25 50

xt 6633 6619 6573 yt 1659 1610 1521 Slip 3177 3246 3397 Ramp angle 43.5 42.3 39.9 4 57.7 56.7 54.6 P/S 2.36 2.35 2.31

spoint 5

10

sdip 2.5

6639 6637 6630

6637 6638 6637

12.5 25 50

1697 1683 1654

1706 1704 1690

12.5 25 50

3155 3159 3178

3186 3181 3150

12.5 25 50

44.0 43.9 43.6

43.5 43.6 44.1

12.5 25 50

58.0 58.0 57.7

57.4 57.5 58.1

12.5 25 50

2.37 2.36 2.36

2.37 2.37 2.36

12.5 25 50

xt 9.0 21.7 21.3 yt 14.9 28.8 26.6 Slip 49.1 64.9 85.0 Ramp angle 0.84 1.0 1.2 4 0.79 1.0 1.2 P/S 0.009 0.015 0.016

5

10

8.1 12.3 20.3

8.4 12.6 17.7

12.0 19.3 30.9

11.9 18.7 28.1

74.2 80.2 98.7

84.7 112 126

1.3 1.4 1.7

1.4 1.9 2.2

1.2 1.3 1.6

1.4 1.8 2.1

0.009 0.013 0.018

0.010 0.015 0.020

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

error in restored dips is not appropriate. Propagating errors will produce a result consistent with the uncertainties in the measurements in the deformed state but with the greater efficiency that inverse modeling achieves over forward modeling. The Waterpocket monocline example shows that our methods are applicable to a natural structure. As with the synthetic section, we found two solutions, one with P/S > 1 and one with P/S < 1. We were able to discard the P/S < 1 solution, however, on the basis of mechanical requirements and analogy to other Laramide folds. This result highlights the fact that while trishear inverse modeling is a powerful mathematical tool, unrealistic solutions are possible. By using a method of inversion that allows us to identify multiple solutions, we are able to discover the full range of kinematicallypossible trishear models, and then to evaluate them on other grounds to narrow the solution further. Our results for the P/S > 1 case are, within reason, similar enough to Cardozo (2005), given the differences between the two approaches. Our standard deviations are relatively small and do not include Cardozo's (2005) model within the estimated uncertainty, despite the fact that this model is visually a reasonable fit to the data (Fig. 12). Thus a strict application of these error estimates may underestimate the range of acceptable models in some cases. The imperfect match between theory (the trishear model), with its simplifications, and the details of a natural structure may be part of the reason for this and is difficult to account for quantitatively. Regardless of the exact uncertainty values that may be appropriate, the results of a test of different data uncertainties (Table 5) indicate that our mean values remain good estimates of the parameter values. Nonetheless, our results for the Waterpocket monocline show advantages to using the APT algorithm as opposed to the grid search and to restoring the data instead of using batch forward models. The APT algorithm allows us with a similar number of models to consider a larger parameter space than the grid search while sampling with much higher precision in the area of interest. For example, when the vertical separation between hanging wall and footwall contacts is well constrained, as in this case and in the synthetic section, there is a clear covariance between total slip and ramp angle. Indeed, our best-fit total slip differs from that of Cardozo (2005) by 627 m, but the uplift of the hanging wall differs by only 11 m due to our correspondingly higher best-fit ramp angle. There is therefore a range of total slip and corresponding ramp angles that will produce approximately the correct uplift, and the grid search finds only the one(s) that happen to lie closest to the grid, while the APT algorithm can search with greater resolution for the best fit and can reveal the covariance between the two parameters. The second benefit of our new approach comes from the dip restoration algorithm we employed, using Eq. (6). The approach used by Cardozo (2005) required the forward modeling of six beds, each of which would consist of a large number of points to be moved individually. If the dips not at formation contacts were to be used as well or if the undeformed elevations were not well known, even more beds would have had to have been modeled, with dips interpolated between them, as in the approach used by Regalla et al. (2010) for the Futaba fault in Japan. By restoring the dips directly, we have only to move and rotate a single point for each data point and we are able to restore equally well dip measurements that do and do not correspond to formation contacts. 5. Conclusions Using the trishear velocity equations of Zehnder and Allmendinger (2000), we have derived an equation for rate of change of dip, Eq. (6), which we have applied to both synthetic and natural datasets. This equation can be used to restore the dips of

171

folded beds by alternately moving the point where the dip was measured and rotating the dip. The error in the restored dip can be used with an estimated uncertainty to calculate probability for each model. Using this approach, dips and points on a bed or contact can be inverted jointly, and each model tested will sample from a probability distribution over the multidimensional parameter space. The traditional grid search algorithm can be used not only to find a best-fit model, but also to estimate the marginal probability densities for each parameter, which can be calculated by integrating over the probability density function. Markov chain Monte Carlo methods, however, provide much faster solution times than the grid search and higher resolution in the areas of interest in the probability density function. In particular, we have found the adaptive parallel tempering algorithm (Miasojedow et al., 2013) well suited to the trishear problem, as it is able to efficiently explore the multi-dimensional parameter space, identify multiple solutions (as represented by multimodal probability distributions), and estimate uncertainty in parameter values. Dip and contact data along a surface profile provide a better constraint on the trishear model than either alone. Dip data, in particular, may provide a poor fit when used alone. It is also necessary to treat errors in dip data in terms of probabilities and to propagate uncertainties through the trishear deformation, as simply minimizing RMS error in restored dips is likely to be inaccurate. Program availability The data inversions described in this paper were performed using a program, InvertTrishear, written for the purpose. The program can be downloaded at www.davidosoakley.com/trishear. html, or can be obtained by contacting the corresponding author. Acknowledgments This research was supported by funding from Geological Society of America Graduate Research grants, American Association of Petroleum Geologists Grants-in-Aid, Sigma Xi Grants-in-Aid of Research (Grant ID #s G20120315159423 and G20130315164439), a Shell Geosciences Energy Research Facilitation Award, and the Penn State Deike grant and Scholten-Williams-Wright Scholarship in Field Geology. We thank Christopher Connors for his helpful review of the manuscript and the editor Joao Hippertt for additional comments. Move software, which was used in the research, was generously provided by Midland Valley under their Academic Software Initiative. References Allmendinger, R.W., 1998. Inverse and forward numerical modeling of trishear fault-propagation folds. Tectonics 17, 640e656. Allmendinger, R.W., Cardozo, N., Fisher, D.M., 2012. Structural Geology Algorithms: Vectors and Tensors. Cambridge University Press, New York. Allmendinger, R.W., Shaw, J.H., 2000. Estimation of fault propagation distance from fold shape: implications for earthquake hazard assessment. Geology 28, 1099e1102. Allmendinger, R.W., Zapata, T.R., Manceda, R., Dzelalija, F., 2004. Trishear kinematic n Basin, Argentina. In: modeling of structures, with examples from the Neuque McClay, K. (Ed.), Thrust Tectonics and Hydrocarbon Systems, vol. 82. AAPG Memoir, pp. 356e371. Bump, A.P., 2003. Reactivation, trishear modeling, and folded basement in Laramide uplifts: implications for the origins of intra-continental faults. GSA Today 13 (3), 4e10. Cardozo, N., Aanonsen, S., 2009. Optimized trishear inverse modeling. J. Struct. Geol. 31, 546e560. Cardozo, N., 2005. Trishear modeling of fold bedding data along a topographic profile. J. Struct. Geol. 27, 495e502. Cardozo, N., Bhalla, K., Zehnder, A.T., Allmendinger, R.W., 2003. Mechanical models of fault propagation folds and comparison to the trishear kinematic model. J. Struct. Geol. 25, 1e18. Cardozo, N., Brandenburg, J.P., 2014. Kinematic modeling of folding above listric

172

D.O.S. Oakley, D.M. Fisher / Journal of Structural Geology 80 (2015) 157e172

propagating thrusts. J. Struct. Geol. 60, 1e12. Cardozo, N., Jackson, C.A.L., Whipp, P.S., 2011. Determining the uniqueness of bestfit trishear models. J. Struct. Geol. 33, 1063e1078. Erslev, E.A., 1991. Trishear fault-propagation folding. Geology 19, 617e620. Erslev, E.A., Koenig, N.V., 2009. Three-dimensional kinematics of Laramide, basement-involved rocky mountain deformation, USA: insights from minor faults and GIS-enhanced structure maps. GSA Mem. 204, 125e150. Hardy, S., Allmendinger, R.W., 2011. Trishear: a review of kinematics, mechanics and applications. In: McClay, K., Shaw, J., Suppe, J. (Eds.), Thrust Fault-related Folding, vol. 94. AAPG Memoir, pp. 95e119. Hardy, S., Ford, M., 1997. Numerical modeling of trishear fault propagation folding. Tectonics 16, 841e854. Hastings, W.K., 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97e109. Judge, P.A., Allmendinger, R.W., 2011. Assessing uncertainties in balanced cross sections. J. Struct. Geol. 33, 458e467. Lin, M.L., Wang, C.P., Chen, W.S., Yang, C.N., Jeng, F.S., 2007. Inference of trishearfaulting processes from deformed pregrowth and growth strata. J. Struct. Geol. 29, 1267e1280. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E., 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087e1092. Miasojedow, B., Moulines, E., Vihola, M., 2013. An adaptive parallel tempering

algorithm. J. Comput. Graph. Stat. 22, 649e664. Pei, Y., Paton, D.A., Knipe, R.J., 2014. Defining a 3-dimensional trishear parameter space to understand the temporal evolution of fault propagation folds. J. Struct. Geol. 66, 284e297. Regalla, C., Fisher, D., Kirby, E., 2010. Timing and magnitude of shortening within the inner fore arc of the Japan trench. J. Geophys. Res. 115, B03411. Smithson, S.B., Brewer, J.A., Kaufman, S., Oliver, J.E., Hurich, C.A., 1979. Structure of the Laramide wind river uplift, Wyoming, from cocorp deep reflection data and from gravity data. J. Geophys. Res. 84, 5955e5972. Stone, D.S., 1993. Basement-involved thrust-generated folds as seismically imaged in the subsurface of the central rocky mountain foreland. GSA Spec. Pap. 280, 271e318. Suppe, J., 1985. Principles of Structural Geology. Prentice-Hall, Englewood Cliffs. Suppe, J., Medwedeff, D.A., 1990. Geometry and kinematics of fault-propagation folding. Eclogae Geol. Helv. 83, 409e454. Tarantola, 2005. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, Philadelphia. Tarantola, A., Valette, B., 1982. Inverse problems ¼ quest for information. J. Geophys. 50, 159e170. Vihola, M., 2012. Robust adaptive metropolis algorithm with coerced acceptance rate. Stat. Comput. 22, 997e1008. Zehnder, A.T., Allmendinger, R.W., 2000. Velocity field for the trishear model. J. Struct. Geol. 22, 1009e1014.