# A geologically-sound representation for evolutionary multi-objective subsurface identification

Vijay Pratap Singh, Marc Schoenauer, Michel Léger

► To cite this version:

Vijay Pratap Singh, Marc Schoenauer, Michel Léger. A geologically-sound representation for evolutionary multi-objective subsurface identification. IEEE Congress on Evolutionary Computation - CEC'05, D. Corne, Sep 2005, Edinburgh - UK, France. pp.454-462. inria-00000855

HAL Id: inria-00000855

https://inria.hal.science/inria-00000855v1

Submitted on 26 Nov 2005

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A geologically-sound representation for evolutionary multi-objective subsurface identification

**Vijay Pratap Singh**
Institut Français du Pétrole
Exploration-Production Business Unit
1 et 4, avenue de Bois-Préau, 92852
Rueil-Malmaison, Cedex France
Email: [V-Pratap.SINGH@ifp.fr](mailto:V-Pratap.SINGH@ifp.fr)

**Marc Schoenauer**
LRI - Bat. 490
Université Paris Sud
91405 Orsay Cedex, FRANCE
Phone: +33 (0) 1 69 15 66 26
Email: [marc@lri.fr](mailto:marc@lri.fr)

**Michel Léger**
Institut Français du Pétrole
Exploration-Production Business Unit
1 et 4, avenue de Bois-Préau, 92852
Rueil-Malmaison Cedex, France
Email: [michel.leger@ifp.fr](mailto:michel.leger@ifp.fr)

## Abstract

A new representation for evolutionary subsurface identification from surface or well geological data is proposed. The idea is to represent the subsurface as the combination of, first, a geologically initial set of horizontal layers and, second, fault parameters like shape and displacements. Based on volume and bed-length preservation, a morphogenesis process then gives the structure at present time. A first implementation of this representation is tested on an artificial geological inverse problem in foothills region: the fault locations and dips are considered as two different objectives, and the ε-MOEA multi-objective evolutionary algorithm is applied. The first results show the efficiency of the chosen representation.

## 1. Introduction

Since petroleum is gradually getting rare, oil companies prospect in new areas like foothills. Because of the rugged topography, the presence of many folds and faults, it is difficult to determine the geologic structures, either by seismic imaging or geological methods.

The problem is hence to identify the geological characteristics of the subsurface from sparse data. The available information can come from either surface geologic surveying or well logging, and consists in layer identification in nature and orientation. The reconstruction of the subsurface thus amounts to an identification problem.

The traditional approach for geological modelling repeatedly evaluates geological models with respect to balancing principles, namely, rock volume and bed-length preservation. These geometrical hypotheses, which are only approximations to reality, mean that

geologic structures behave locally like a folded paperback book, with layers slipping on each other like pages during folding. This geological modelling, which may be backward (Moretti et al., 1989) or forward (Endignoux et al., 1989), is very cumbersome and human time consuming.

Geophysical modelling, on the other hand, has taken advantage from the progresses of numerical modelling, and several works address the problem using classical gradient-based identification methods. Because the problem either has to be oversimplified, or is ill-posed, several works used Evolutionary Algorithms (EAs), that are known to be well suited for solving complex, mathematically ill-posed identification problems (Schoenauer and Sebag, 2002). However, a critical question is the choice of a representation (and of the associated variation operators (crossover and mutation)). In the case of subsurface identification, state-of-the-art works either assumed some expert knowledge about the geometry of the subsurface structure (e.g. horizontal layers, or more complex geologically-driven topology (Boschetti et al., 1996)), or used global automatic models that lack geological soundness (Mansanné et al., 2002) (e.g. the subsurface structures are not balanced).

In this paper, we propose an original indirect representation for subsurface identification. The morphogenesis process (from the encoded representation to the present subsurface structure) only generates balanced geometries. Moreover, it can be used to solve either or both the geological and the geophysical identification problems by EAs automatically, and without the need for any expert
knowledge. The idea is to represent a subsurface
structure by the combination of, first, an initial
configuration (series of flat homogeneous layers of
different geological nature), second, some initial faults
in that simple structure, and third, a description of
successive geological deformations along the ages. A
kinematic model, that relates fold geometry to fault
shape and displacement, is used to construct cross-
sections of the geological structures. These cross-
sections are automatically balanced by virtue of the
above-mentioned balancing principles (volume and
bed-length preservation) upon which the kinematic
model is based. The resulting complex subsurface
structure is then used to compute the identification
criteria: comparison with surface or subsurface
geological data, geophysical simulation of some
seismic experiment, etc.

However, as already mentioned, the different types
of available information may be conflicting. Therefore
the problem is better reformulated into a multiobjective
optimization problem: there is no unique optimal
solution but a set known as the Pareto set of optimal
solutions, also termed non-dominated non-inferior,
admissible, or efficient solutions. In this paper,
multiobjective evolutionary algorithms (MOEAs),
more precisely, the Epsilon-MOEA [Deb-epsilon],
have been used to tackle this multiobjective problem.

The paper is organized as follows: first the
geological framework is presented in section 2,
together with the kinematic model and the complete
geological inverse problem is described. Section 3
introduces the original indirect representation and the
variation operators, as well as the different fitness
criteria. Section 4 gives the first experimental results
on two artificial geological problems, involving
respectively 3 and 7 faults, and discusses the limitation
of such purely geological experiments. Finally, section
5 sketches the directions of research that the proposed
approach opens.

## 2. The geological problem

2.1 Geology

Complexity and ambiguity of mountain fronts are
what makes them the most challenging fronts for a
geologist or a geophysicist. Now the thrust regions
(mountainous regions) are the next human quest for
petroleum exploration. These fold-and-thrust belts are
the manifestation of collision of tectonic plates. Most
fold-and-thrust belts consist of foreland-progressing
and hinterland-dipping thrusts (Boyer and Elliott,
1982) that may connect layer-parallel detachments
(décollement horizons). The distribution of mesoscopic
deformation features in fault-related folds has
important consequences for hydrocarbon migration,
trapping and production. Structural, environmental and
stratigraphic factors control the spatial distribution and
intensity of the deformation features (e.g. Fischer and
Jackson, 1999) of mountains. Under a given
rheological stratigraphy (Woodward and Rutherford,
1989) and prevalent environmental conditions (Stewart
and Alvarez, 1991; Jamison, 1992; Lemiszki et al.,
1994). It is fault-fold kinematics which controls the
distribution and intensity of deformation that develop
in a rock sequence (Fischer et al., 1992; Storti and
Salvini, 1996).

**2.2 Kinematic Models**

For the representation of mountainous structure
mainly three types of kinematics models appear in the
literature. 1) kink model (Suppe 1983, 1990), 2) trishear model (Erslev, 1991; Hardy and Ford, 1997; Allmendinger, 1998) and 3) force fold model (Johnson and Johnson, 2000). Amongst the above mentioned three types the kink model is considered best suited for thin skinned tectonic regions. Since we are interested in such regions, hence kink model is being used.

A two-dimensional geometric model of fault-bend
folding for a thin-skinned tectonic was first formulated
by Suppe (1983) based on the conservation of layer
thickness and bed length. He showed its applicability in
the simulation of the Pine Mountain thrust sheet in the
southern Appalachians, and to the fold and thrust belt
of Western Taiwan. The use of kink-band style folding
with non-deforming footwalls (characteristic of such
geometric models) has been criticized by Ramsay
(1991) as being unrealistic. Another fault-bend folding
model was proposed by Contreras and Suiter (1990),
however it was also based on the conservation of area.
Since their introduction, these geometric models have
received a lot of attention (Medwedeff, 1989; Jordan et
al., 1993) and have been used extensively to predict
thrust fault geometries at depths based on observed
fold geometries. They have also been used in a forward
modeling sense to predict hanging wall geometries
above thrust faults (Mitra, 1990; Contreras and Suiter,
1990, 1997; Zoetemeijer, 1993; Hardy S., 1995).

## 2.3 The Contreras Model

The kinematic model used in the current simulation
was introduced by Contreras and Suiter (1990). This
model can be applied to regions being deformed by
shortening, extension and also for the duplex systems
(Contreras and Suiter, 1997). This is a two-dimensional
model and is based on a coordinate transformation
from less deformed state to a more deformed state. The
medium is subdivided into domains of constant dip and
homogeneous displacement vector fields that are
delimited by the planes bisecting the fault inflections.
The displacement occurs by translation. The
displacement trajectory is of constant length for all the
displaced particles throughout the medium and is
parallel to the underlying active fault segment. This
model also considers fault parallel simple shearing. A
complete derivation Contreras fault bend folding (FBF)
was provided by Contreras and Suiter (1990). The
heterogeneity of the displacement vector field across
axial planes introduces longitudinal and angular shear
strains. Transformations from one state to another state
do not cause a change in area and also the deformation
is isochoric, a characteristic typical of deformation by
simple shear (Truesdell and Toupin, 1960).

# 3. The Evolutionary Algorithm

## 3.1 The Representation

As mentioned in the introduction, the main difficulty
when tackling a geological or geophysical
identification of subsurface structure problem by EAs
is to define a representation: structured representations
are either unrealistic in complex regions like foothills
(e.g. a fixed number of horizontal layers) or require
some very specific expert knowledge about the
geometry of the subsurface structure (Boschetti et al.,
1996). Unstructured representations, on the other hand,
don't require any input from the expert – but they lack
of the minimal geological common-sense and hence
can lead to quasi-perfect numerical fit that are absurd
from a geological point of view (Mansanné et al.,
2002). Note that some alternative possibilities have
been proposed, that define a global model with
numerous local parameters (e.g. using spline nodes
uniformly spread on the considered domain), but such
representations have the same lack of geological
soundness, while additionally resulting in a huge
optimization problem for the sake of precision.

Moreover, it seems very difficult indeed to try to
constrain the solution proposed by unstructured
representations with some geological rules (e.g. some
simple rules could be that the underground velocity
should increase with the depth – but interesting regions
are precisely regions where such rule is violated!).

The representation proposed in this work deals
with the above-mentioned problem by relying on a
kinematic model of subsurface deformation from an
initial simple state (horizontal layers, that can be
assumed to be universally true in all regions at some
point in the past) and subsequent deformations based
on one of the kinematic models discussed in section 2.2
that ensure the consistency of the structure at any time,
including the final state of the model – the target
subsurface structure.

In EAs terminology, the **genotype space**, in which
the algorithm will actually search, and where the
variation operators are defined, is here the space of
initial configurations plus initial faults plus series of
deformations (see figure 1 for the detailed description
of the complex data-type). The **phenotype space**, or
behavioral space, where the fitness of each genotype is
computed, is the space of subsurface structures,
obtained from the genotypes by applying a
transformation, called the **morphogenesis process**, or
also the numerical Contreras model described in
section 2. The tectonic modelling presented in section 2
is sound because we will obtain only balanced
subsurface structures by this morphogenesis process.

<figure><img src="image_2.png" /><figcaption>**Figure 1** *A chromosome: Unknown parameters are the number of initial layers and, for each layer, its thickness; the number of faults and, for each fault, its origin, and number of segment, and, for each segment, its length and dip; and the displacements (for each fault, one horizontal displacement).*</figcaption></figure>

**3.2 Initialization and Variation Operators**

First note that the representation proposed in this work
is possibly a variable length representation (variable
dimension parameter space). However, in the first
experiments presented in section 3, we have used a
fixed length representation: the number of faults and
the number of segments per fault will be fixed (from
the target model). We shall hence here only describe
the operators of this fixed-length restriction.

Initialization

Here a complex structure of real values is being used,
therefore the parameters are initialized uniformly on
a given interval, with some geometrical constraints:

* Displacement of a fault can not be more than the difference of distance between the first segments of two successive faults.

* The length of first segment increases in the successive faults.

Crossover Operator

Because the representation is at the moment fixed
length, the Simulated Binary Crossover (SBX) (Deb
and Agrawal 95) recommended with the epsilon-
MOEA algorithm has been used. During the crossover,
the geometry constraints listed above are respected.
A quick parameter study (12 independent runs for each
value of η) has been performed (on the 3-faults
problem described in section 4). All final populations
are merged and non-domination sorting is performed
(result can be seen in Figure 2). It is found that high
value of η is most suitable for this simulation, and
thereon the distribution index η has been set to 15.
Note that the recommended values of ε₁= 0.05 and ε₂=
0.05 were found robust enough and used in all
experiments presented in section 4.

<figure><img src="image_3.png" /><figcaption>**Figure 2** *Influence of the distribution index η of SBX crossover on the results (results of 12 runs).*</figcaption></figure>

However, because the SBX crossover can hardly be
applied on variable-length genotypes, some more
variable-length specific crossover operators are used in
the variable-length case (not described here for space
reasons).

Mutation Operator

The SBX crossover is known to be self-
adaptive in the sense that the spread of the possible
offspring solutions depends on the distance between
the parents, and decreases as the population converges
(Beyer and Deb 2001). Hence it is generally used
without any mutation (the use of mutation together
with SBX in the context of the problems presented in
section 4 did not seem to bring any benefit indeed, and was abandoned).

However, because the variable-length case requires additional use of mutation, self-adaptive Gaussian mutation (à la ES) is used on all real-valued parameters, while specific variable-length mutation operators are used to modify the numbers of faults and segments in each fault.

### 3.3 The ε-Multiobjective Evolutionary Algorithm

As mentioned in the introduction, various sources of data will be used for identification. Hence, the identification problems will be turned into Multi-Objective Problems (MOPs). Examples of such problems can be seen in section 4.2 for the geological modeling problem, or in [Mansanné] for purely geophysical identification using both the semblance and the least-squares error on seismic signals. But the ultimate goal of the present work is to simultaneously use all available data, both geological and geophysical, leading to even more objectives.

During the past decade, various multiobjective evolutionary algorithms (MOEAs) have been proposed. A comprehensive review are given by Coello et al., (1999); Veldhuizen (2000); Deb (2002). Although these MOEAs differ from each other in both exploitation and exploration, they share the common purpose: searching for a near-optimal, well-extended and uniformly diversified Pareto optimal front for a given MOP. For the optimization of foothill structure (section 4), we have chosen the recently proposed, ε-MOEA (Deb et al., (2003)) because it achieves all the above desired properties at a minimum computation cost. The code of this algorithms is graciously offered by the authors at URL http://www.iitk.ac.in/kangal/soft.htm.

ε-MOEA is based on the principle of ε-dominance, that relaxes the strict Pareto-dominance principle. It uses two co-evolving populations: an EA population and an archive population A, and is a steady-state EA. The initial archive population A is assigned from the initial EA population using ε-non-dominated solution.

Thereafter two solutions are selected, one from the EA population using tournament selection and the other from archive A using random selection. From those parents, an offspring O is created, and is used to update the EA population on the basis of strict dominance, and the archive population on the basis of ε-dominance. Only one solution in each hyper-box on the Pareto-optimal front is allowed, to favor diversity.

## 4. First Results in Geological Modeling

### 4.1 The Geological Identification Problem

For identifying a foothill structure, a purely geological inverse problem is defined. In this inverse problem, the unknowns are the parameter of a structure, and the data are layer dip or fault location measured on the topographical surface or along a well (Figure 3). Experimental results are being presented on artificial problems: the “experimental” data are computed on a model that is in itself the result of an evolution using the Contreras kinematic model starting from a known initial configuration. Two models have been considered so far, a 3 fault and a 7 fault models.

<figure><img src="image_4.png" /><figcaption>Figure 3 Data are the measured fault location (dark circles) and the dip of the layers (light circle) along surface and well.</figcaption></figure>

### 4.2 The Evaluation Functions

This is the most costly and important step for real applications. Contreras and Suiter (1990), kink model is used for the deformation of the model from the initial state to final state. After deformation, dip and fault location parameters are evaluated for optimization purpose. In general, dip of faults, stratigraphic unit, and fault location are easily obtained from the field observations, well logs (dip-meter) or by remotes sensing data. For the evaluation of a model, a least-squares criterion measures the discrepancy between the field data and corresponding dip and faults location on
the current model. The relation is shown by the following equation.

$$f1 = \sqrt{\sum_{i=0}^{nd} (Dip\_ob_i - Dip\_gen_i)^2} \\ f2 = \sqrt{\sum_{i=0}^{nf} (Fault\_position\_Ob_i - Fault\_position\_gen_i)^2}$$

The computational cost of the evaluation function is hence negligible when compared to that of the morphogenesis process (the Contreras model that computes the deformation of the subsurface structure). For instance, for the experiments presented in next subsections, the total computational cost on a Pentium 3.4GHz is about 4 hours.

## 4.3 A Three Fault Model

A target mountain front model, with five layers and three faults, is shown in figure 4. Because this model was generated from an initial genotype (as described in section 3.1) with the Contreras model, it is the exact solution to the identification problem. On this example, data is only available at the surface, light circles representing dip information, dark circles fault locations.

<figure><img src="image_5.png" /><figcaption>**Figure 4** Good results on the 3-faults model: (a) is the target structure. (b) and (c) are result structures. Structure (b) is similar to target structure in all respects. Structure (c) is also similar to the target structure but fault locations are slightly different. However, both structures coincide with the target at the observation points.</figcaption></figure>

For this experiment, the population size was 30, archive size 30, and number of generations 300, and, as mentioned in section 3, the representation was restricted to fixed-length 3 fault genotypes.

Figure 4 shows some identified models that are very similar to the target model. Model 4(b) is similar to the target model in terms of surface and subsurface fault location and topography, while model 4(c) is also similar to the target structure but its fault positions are slightly different.

However, some others models have also be obtained, which are numerically quasi-optimal too, but visually quite different from the target model (see figure 5). This leads to the conclusion that the problem is underdetermined: we need more information in order to uniquely define a solution.

<figure><img src="image_6.png" /><figcaption>**Figure 5** (a) and (b) are the others optimized models. The models (a) and (b) are similar along the topographic surface, but visually different below.</figcaption></figure>

## 4.4 Results on the Seven Fault Model

A more complex seven fault model is shown in figure 6(a). On this model the observation points are located both on the surface and in a well. For this experiment population size was 50, archive size 50, and number of iterations was 500. The identified models are shown in figure 6(b) and (c). Again, they are very similar to the target model. In figure 6(b), the fourth fault is visually different from its analogue in the target model but the rest of the faults are both numerically and visually similar. Figure 6(c) also looks very similar to the target but for the sixth and fourth faults. Here also, as can be seen on figure 7, some identified models are quite
different from the target model – though being quasi-optimal numerically, i.e. the available data on the surface and in the well are almost perfectly fitted.

This situation is very typical of ill-posed problems: we do not have sufficient information to ensure the uniqueness of the solution.

Remember that some basic geometry constraints have been added to the representation (e.g. successive faults first segment length will be more than the previous one, see Section 3.1). Such additional knowledge did prove beneficial: without this information, the algorithm need around 1000 iterations to converge while only around 300 with the constraints. However, higher order of information like seismic, gravity or magnetic data, seem necessary to improve the results.

<figure><img src="image_7.png" /><figcaption>**Figure 6 (a) is the target structure, (b) and (c) are the numerically optimised structures. In spite of minor differences, models (b) and (c) are fairly similar.**</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7 (a) and (b) are other numerically optimal structures for the 7-fault model. But models 7(a) and 7(b) are obviously different from the target model 6(a). However they fit rather nicely available dip and fault position data as well.**</figcaption></figure>

# 5. Conclusions and further directions

The complexity and ambiguity of mountain fronts pose a significant challenge for both geologist and geophysicist in determining the geological and geophysical parameters. In traditional approach mountain front identification needed either rigorous human interaction because of absence of established numerical criteria to act as an inversion target, noisy surface and subsurface data and little prior information, or lead to geologically unrealistic results because, except for very simple regions, the geometry of the subsurface structure cannot be accurately predicted.

We have proposed a new representation for subsurface structure that has both advantages: it is geologically relevant (all structures are balanced), and it does not require human guesses about the geometry of the unknown structure. Moreover, this is the first time, to the best of our knowledge, that multiobjective optimization is applied to a subsurface identification problem. The first results, obtained on synthetic geological identification problems, show the power of the proposed representation, even if we restricted this representation to fixed length for those preliminary experiments.

However, those first results also made it obvious that the geological identification problem is ill-posed, and that more data is needed in order to reach good solutions with more robustness. On-going work is indeed dedicated to adding the seismic data that are available. However, the computation of the geophysical part of the fitness will now rely on some heavy numerical simulation and the overall computational cost will be greater by one order of magnitude. One critical issue will then be how to balance the workload between the purely geological and the geophysical objectives.

Another direction is of course to experiment the full power of the proposed representation with respect to the variable number of faults: at the moment, though we claim that no expert knowledge is needed, we did
specify the correct number of faults because we used synthetic data. Though there are some good hints about the number or faults thanks to well information, there is always a chance that we miss some faults, or fail to connect the same faults in two different wells – and using the variable-length representation will then make full sense.

**References:**

1. Allmendinger, R.W. (1998) "Inverse and forward numerical modeling of trishear fault-propagation folds". *Tectonics* 17 (1998), 640–656.

2. Beyer Hans-Georg, Deb K (2001) "On self-adaptive features in real-parameter evolutionary algorithms" IEEE Trans. Evolutionary Computation 5(3): 250-270

3. Boyer, S.E., and Elliott, D. (1982) "Thrust systems": American Association of Petroleum Geologists Bulletin, v. 66, 1196-1230.

4. Boschetti, F., Dentith, M.C. & List, R.D. (1996) "Inversion of seismic refraction data using genetic algorithms", *Geophysics*, 61, 1715–1727.

5. Chester J.S. (2003) "Mechanical stratigraphy and fault-fold interaction, Absaroka thrust sheet, Salt River Range" Wyoming, *Journal of Structural Geology* 25 (2003), 1171–1192.

6. Contreras, J., Suiter, M. (1990) "Kinematic modeling of cross-sectional deformation sequences by computer simulation". Journal of Geophysical Research 95, 21913-21929.

7. Contreras J. (1991) "Kinematic modeling of cross-sectional deformation sequences by computer simulation: coding and implementation of the algorithm", Computers & Geosciences, v.17 n.9, 1197-1217.

8. Contreras, J., Suter, M. (1997) "A kinematic model for the formation of duplex systems with a perfectly planar roof thrust". Journal of Structural Geology 19, 269-278.

9. Coello C.C.A. (1999) "A comprehensive survey of evolutionary-based multi-objective optimization techniques" *Knowledge and Information Systems An International Journal* 1 (3), 269–308.

10. Dahlstrom C.D.A. (1990) "Geometric constraints derived from the law of conservation of volume and applied to evolutionary models for detachment folding" *American Association of Petroleum Geologists Bulletin* 74 (1990), pp. 336–344.

11. Deb K. and Agrawal R. (1995). "Simulated binary crossover for continuous search space" *Complex Systems*, 9(2):115-148.

12. Deb K. (2001) "Multi-Objective Optimization using Evolutionary Algorithms" Wiley, Chichester, UK.

13. Deb, K, Mohan, M. and Mishra, S. (2003). "A Fast Multi-objective Evolutionary Algorithm for Finding Well-Spread Pareto-Optimal Solutions" *KanGAL Report No. 2003002*.

14. Endignoux L, Moretti I., and Roure F., 1989, "Forward modeling of the southern Apennines", Tectonics V8, 5, 1095-1104.

15. Erslev E.A. (1991) "Trishear fault-propagation folding" *Geology* 19 pp. 617–620.

16. Fischer, M.P., Jackson, P.B., (1999) "Stratigraphic controls on deformation patterns in fault-related folds: a detachment fold example from the Sierra Madre Oriental, northeast Mexico" Journal of Structural Geology 21, 613–633.

17. Hardy, S. (1995) "A method for quantifying the kinematics of fault-bend folding" *Journal of Structural Geology* 17, 1785-1788.

18. Jamison W.R. (1992) "Stress controls of fold thrust style" *Thrust Tectonics*, Chapman and Hall, London (1992), 155–164.

19. Jamison W.R. (1997) "Quantitative evaluation of fractures on Monkshood anticline, a detachment fold in the Foothills of Western Canada" *The American Association of Petroleum Geologist Bulletin* 81 (1997), 1110–1132.

20. Mansanné, F., and Schoenauer, M., (2002), "An automatic geophysical inversion procedure using a genetic algorithm" Soft computing for reservoir characterization and modeling : Physica-Verlag, Heidelberg, 331–335

21. Mitra, S. (1988) "Three-dimensional geometry and kinematic evolution of the Pine Mountain thrust system, Southern Appalachians" Geological Society of America Bulletin, v. 100, 72-95.

22. Mitra, S. (1992) "Balanced structural interpretations in fold and thrust belts" Structural geology of fold and thrust belts, Johns Hopkins University Press, Baltimore, p. 53-77.

23. Moretti I. and Larrère M., 1989, "LOCACE: Computer-Aided Construction of Balanced Geological Cross-sections", Geobyte, Oct. 1-24.21.

24. Morgan J.K.and Karig D.E. (1995) "Kinematics and balanced cross-section across the toe of the eastern Nankai accretionary prism" *Journal of Structural Geology* 17, 31–45.

25. Ramsay, J.G. (1992) "Some geometric problems of ramp-flat thrust models" Thrust tectonics, Chapman-Hall, London, 191-200.

26. Rich, J. L. (1934) "Mechanics of low-angle overthrust faulting as illustrated by Cumberland thrust block, Virginia, Kentucky and Tennessee" Bulletin of the American Association of Petroleum Geologists 18, 1584-1596.

27. Suppe, J. (1983) "Geometry and kinematics of fault-bend folding" American Journal of Science, v. 283, p. 684-721.

28. Suppe, J. and Medwedeff, D.A. (1984) "Fault-propagation folding" Geological Society of America Abstracts with Programs, v. 16, 670.

29. Suppe, J. and Medwedeff, D.A. (1990), "Geometry and kinematics of fault-propagation folding" Eclogae geologicae Helvetiae, v. 83, 409-454.

30. Tanner, P.W.G. (1992) "Morphology and geometry of duplexes formed during flexural-slip folding" Journal of Structural Geology, v. 14, 1173-1192.

31. Sen, M. K. and P. L. Stoffa, Nonlinear one-dimensional seismic waveform inversion using simulated annealing, *Geophysics*, 56, 1624-1638, 1991.

32. Schoenauer M.and M. Sebag M, (2002), "Using Domain Knowledge in Evolutionary System Identification" Evolutionary Algorithms in engineering and Computer Science John Wiley

33. Storti F. and Salvini F. (1996) "Progressive rollover fault-propagation folding: a possible kinematic mechanism to generate regional-scale recumbent folds in shallow foreland belts" *American Association of Petroleum Geologists Bulletin* 80, 174–193.

34. Veldhuizen D. and Lamont G. (2000) "Multiobjective evolutionary algorithms: analysing the state-of-the-art" *Evolutionary Computation* 8 (2), 125–147.

35. Woodward N.B. (1999) "Competitive macroscopic deformation processes", *Journal of Structural Geology* 21, 1209–1218.

36. Zoetemeijer, R. (1993) "Tectonic Modeling of Foreland Basins: thin skinned thrusting, syntectonic sedimentation and lithospheric flexure" Ph.D. thesis, University of Amsterdam.