GeoComputation Logo

GeoComputation 2000 HomeConference ProgrammeAlphabetical List of Authors

Computer simulation of hangingwall deformation

Sergio Perez
Geology Department, Royal Holloway,University of London, Egham Hill, Egham, TW20 OEX, Surrey, U.K.

1. Introduction
2. Programming the model
2.1. Set up of initial conditions
2.2. The model is posed and solved
2.2.1 Analysis of the derivatives of the inclined simple shear function
2.3. Interpolation of  the velocities
2.4. Graphical display of the simulation


A computer simulation based upon a finite difference model of hanging wall deformation is elaborated via Mathematica Programming.

Vertical simple shear is usually adapted to set up the boundary conditions because it facilitates the calculations.

Nevertheless, in the program is used  a Mathematica®-built displacement direction function consistent with inclined simple shear conditions.

To find the derivatives of the displacement direction function, an essential sub-routine in the simulation, two approaches are used. The first involves a numerical solution based upon the Crank-Nicholson method. The second requires the calculus of the derivatives with the support of Mathematica® built-in functions.

A brief analysis is carried out on  the graphical display of these sub-routines.

An assesment on the efficiency of Mathematica® for linear finite-element modelling is given based upon the availability and efficiency of resources found in the system during the programming  of the algorithms.

The graphical and analytical capabilities of Mathematica® are used to generate the trajectories and velocity fields of the hanging wall particles with respect to an Eulerian system of reference.

Key words: Hanging Wall Deformation, Simple Shear, Finite Difference, Crank-Nicholson.

1. Introduction

Forward modelling of hanging wall deformation can be based upon finite-difference solution of a partial differential equation (Waltham, 1989).

Volume balance is achieved as a consequence of constraining the deformation to incompressible conditions, but the effects of compaction (Waltham, 1990) and density increase with depth (Waltham, 1992) can be included to support more elaborated models.

This forward modelling is specially useful for the interpretation of seismic data from fault zones, that are known to be of poor resolution (Kerr and White, 1994).
Comparing the simulation of deformation associated with a variety  of fault shapes with either  the deformation of the overburden above the fault zone  or the one of the faulted sequence  makes feasible to deduce the real fault geometry, partial o totally faded in the seismic record.

For oil business purposes this branch of computer modelling is now of increase demand since makes possible valid and reliable interpretations of the boundaries of faulted reservoirs (Kerr and White, 1996). The impact of this fact is felt specially for a reliable calculation of the volume of such reservoirs and  for the plans of well drilling toward its boundaries.

2. Programming the model

The basic ideas of the method designed by Waltham (1989) and the Mathematica®-based simulation are:

2.1. Set up of initial conditions

One of the first tasks of the simulation is to pre-define a fault over which hanging wall deformation would occur. The dip of the fault at the point of intersection with a pre-determined isogon provides the displacement direction of the velocity for any point along the line (Figs. 1.1 and 1. 2).

 The coordinates of the nodal points in the hanging wall can be arranged in tables using the respective Mathematica® facilities.

  The points along the isogons (-20 degrees of  slope) would have a displacement direction parallel to the tangent to the fault at its intersection. In Figure 1.1 are represented the points along the isogon y=5.67(x-4)+2 (black dots) and the isogon at the boundary between the horizontal line (detachment) and the curved section of the fault . That isogon correspond to the line y=5.67(x-10). The line separates the  displacement directions of zero angle  from those continously decreasing along the circumference. The fault would sole out in a horizontal plane above which horizontal layering would be  represented joining the plot of a nodal set of points prior to deformation (Fig. 1.2). The horizontal line between the ordinates 3 and 4 would be a detachment surface during inversion (see Fig. 7.2)

Mathematica® allows to obtain the displacement direction for any point in the hanging wall by the means of a piecewise function, defined analytically following the conditions of inclined simple shear .

That function is subsequently used for the calculations involved in the Crank-Nicholson method and to define, later on, the position of hanging wall nodes as its deformation evolves.

From Figure 2.1. is possible to identify a zone of subtle torsion in the surface. Is located  around the boundary between the horizontal line and the circumference. In the model  corresponds to the change between inclined simple shear angle domains, from 0 radians  to the region where depends upon the slope of the tangent to the circumference.

The real boundary between the displacement directions domains is an inclined simple shear line. In Figure 2 is the line at the top of the vertical plane, outlined in red,  above the surface. That line corresponds to the equation
y =5.67(x-10), and is the inclined simple shear line at the intersection between the detachment and the listric section of the fault.

The region of  torsion in the surface is in reality an artifact due to the  combination of two factors:

1) The inclination of the simple shear lines does not make the correspondence of angular values for the points in the domain parallel to the Y-axes.

2)  For graphical purposes, Mathematica® uses a default grid size, and interpolates the function between  the points located in the grid. As a matter of fact, the Mathematica® default sample interval cross the real boundary (y=5.67(x-10)) between the displacement direction domains, as is seen in Figure 2.1.

A better graphical approach to the function is reached if the plot is made based upon a smaller interval size for the grid than the default one. In Figure 2.2. is appreciated that the smaller interval grid size allows an approach of the plot to the real surface boundary at y=5.67(x-10) that is closer than the one given by the default grid.

2.2. The model is posed and solved

The  solution of the system of lineal equations with the unknown velocities of nodal points  is obtained through the following sequence of calculations:

The boundary conditions for the velocity field are set up following Walthamís guidelines. Of particular interest is the fact that at the boundary lateral side of the model a constant horizontal velocity, either of extension or normal inversion, is endowed.

The magnitude of the unknown  velocity field is calculated from the system of linear equations generated from the partial differential equation that governs the model. In other words that implies to find a solution of the Laplace Equation for potential flow (Harbaugh and Bonham-Carter, 1970; Waltham, 1989):

2.2.1 Analysis of the derivatives of the inclined simple shear function

One of the feasible solutions of the Laplace Equation follows the finite-difference scheme provided by Waltham (1989).

To that end the Mathematica® built-in function for derivatives allow a system-incorporated  mean to get the derivatives of the displacement direction function illustrated in Figure 2, either along the X axis or along the Yaxis (Figs. 3.B and 3.D, respectively). To codify the Crank-Nicholson solution (Kreyszig, 1995), as suggested by Waltham (1989), is  also a feasible task (Figs. 3.A and 3.C).

The difference between the graph of the functions is seen at the interval where there is a torsion on the displacement direction surface, as is depicted in Figures 3.A to 3.D . The Crank-Nicholson scheme produce a smooth torsion surface at the transition while the built-in function produce a sharp one for the same default interval grid size, common to all  plots.

Again, the default grid size that the system uses to make plots should be refined to improve the insight in the real behavior of the functions. In Figures 4.A to 4.D the plots of the derivatives obtained  by the Crank-Nicholson scheme give a steeper surface at the zone of the boundary between the displacement angle domains.

The plots for each direction of the derivatives look more similar than what is appreciated in Figure 3. Therefore, the difference between the plots seen in Figure 3 is consequence of the grid used. As a matter of fact, as the  functions approaches asymptotically the real value of the derivative,  the plot of both derivatives along each direction should coincide as the grid is refined.

2.3. Interpolation of  the velocities

 As the finite-difference solution of the Laplace Equation provides the  rate of displacement for a particle at any point of the hanging wall, its deformation is finally modelled in combination with the known displacement direction at  successive intervals of time.

This process involves the calculation of velocities for points outside the grid. This is  consequence of the displacement of the nodal points  as deformation evolves.
 A linear interpolation function is suggested to be used to find the velocities of the displaced points.

To tackle this is necessary to find a way to make interpolation over an scatter, non-regularly spaced distribution of points.

To the authorís knowledge only the ExtendedGraphics package developed by Wickam-Jones (1994) provides an answer to this situation. The TriangularInterpolate function from this package does allow to make linear interpolation of scattered distribution of points.   It can be downloaded  from the  Mathematica® website.

The graph of the function applied to the extension and inversion velocity data is seen in Figures 5.1 and 5.2, respectively. The Figure 5.3 displays the default delaunay triangulation, in plan view, that  is (inferred) to be used by the TriangularInterpolate function.

As it is appreciated in Figure 5.2, the plot of the dataset obtained from the TriangularInterpolate function (made using  the built-in function ListSurfacePlot3D) assumes a triangulation scheme that interpolates points from the listric fault across the upper ramp-detachment level, and leaves an open, non-interpolated narrow area between the ordinates 4 and 3.

But one thing seems to be the triangulation followed by the display of the function (Figure 5.2) an another the triangulation followed by the TriangularInterpolate(Figure 5.3) function to make the interpolation. Since both are not necessarily the same, one way to tackle this inconvenient  is to split the  interpolation  in two zones where the TriangularInterpolate function and ListSurfacePlot3D  operate according to the real boundary domain of the model (Figs. 6.1, 6.2 and 6.3).

At this stage, the rate and direction of displacement for a particle at any point of the hanging wall during extension or inversion can be displayed using the function for visualising vectorial fields (Figs. 7.1 and 7.2):

2.4. Graphical display of the simulation

 The desired computer model is obtained using a convenient time-based loop including the visual display of the deformation  of  the bedded sequence (Click here to see an animated sequence of the simulation) .
 As Waltham (1989) pointed out, the model allows to get insights on listric faulting evolution and syntectonic sedimentation during extension (first frames of the animation).

Nevertheless,  the basic algorithm programmed does not include neither the difussion equation nor the compaction due to burial (Waltham, 1992). Hence is an ideal simulation of syn-tectonic sedimentation below water level at  shallow burial depths.

During inversion is feasible to analyse the dynamics of features like the evolution of the position  of the null point, the development of the  inversion anticline coupled with the fault-related fold  and the deformation of the pre and syn-inversion sequences (last frames of the animation).


The results presented in this paper introduces Mathematica® programming for a basic finite-difference model of hanging wall deformation that actually accounts with more elaborated treatment and algorithms programmed in Fortran and  C languages.  Nevertheless the high level of technical resources supplied by Mathematica® in its own computing environment and its capacity to link with Fortran and C packages implies a short term leap from fully developed programs in those languages to Mathematica® programming.

One of the features that the finite-difference model designed by D. Waltham awaits for enhancement deals with 3D out-of-plane movements and deformation of the hanging wall, like happens during transpression or transtension  in strike-slip settings.
The resources provided by high level systems of doing mathematics by computing seems promising to supply  tools that could pave over the way toward developments in such topics.

Respect to Mathematica® itself,  two comments are made.

First, the brief inspection of the outcomes for the derivatives of the angle of displacement direction indicates that care must be taken to crosscheck the effects on the calculations due to the default interval sample that the system uses when operates over functions.
The desired level of accuracy can be obtained changing the default settings accordingly.

Second, the system is under active development by a Research Group and hopefully high standard built-in  tools will be available soon to tackle efficiently all the spectra of specific problems that arise in the process of interpolation of non-uniformly datasets, usually found to develop  finite element  models, either linear or non-linear.

But by now the ExtendedGraphics package is unsupported and seems to have inconvenient handling figures with the default format from lists, so undetermined values could be returned from points within the interpolated area during the run of the loops.  Step-by-step calculations are then needed to make the interpolations and the process is time consuming.

The default triangulation process of the TriangularInterpolate function is difficult to visualice  for complex shapes usually found at the boundaries of geological models, but not to say finite element domains in general. One standing question is whether or not ListSurfacePlot3D makes the same triangulation as the TriangularInterpolate function. Splitting the data to create  piecewise interpolating functions is a solution to avoid doubtful results, but hopefully a more advanced version of  ListSurfacePlot3D will be available soon to control visually the boundaries where  the TriangularInterpolate function  or an equivalent one is needed to be adjusted.


1. Mathematica® as an efficient  derivative modelling  tool and its capabilities for symbolic algebra and calculus are feasible to use for computer simulation of algorithms based upon linear finite-difference designs as the one described by Waltham (1989). The development of the  linear finite-element program is enhanced by the use of the built-in functions and operations on lists supplied by the system.

2. Regarding the resources required from the system  for linear finite-element programming, a wise and critical attitude should be maintained over  what is  Mathematica® doing. Specifically to  avoid a "black box" usage of   the tools supplied by Mathematica®  it is recommended that:

2.1. When  the graphical outcome of  functions are analysed, care must be taken to distinguish between what is consequence of the mathematical properties of the function itself and what is consequence of interpolation over the grid interval size applied to its graphical display.

2.2.When complex boundary shapes are involved in the modelling it is recommended to crosscheck the default triangulation process that operates for ListSurfacePlot3D if this function is applied to visualise values given by the function TriangularInterpolate. The default triangulation process that operates for  TriangularInterpolate itself should be controlled as well.
The application of those built-in  functions over domains splitted off non-complex boundaries would  be a confident way  to use them.


This project is supported by a post-graduate scholarship the author receives from PDVSA, Gerencia de Tecnologia.

The author wish to give thanks to the following persons:

To Dr. Dave Waltham, by his valuable comments and suggestions.

To Cipriano Cruz, now at the Coordinacion Academica of the U.C.V., by the benefits the author owe to his standard of excellence in undergraduate calculus teaching.