Dynamic restoration of profiles across diapiric salt structures: numerical approach and its applications

Tectonophysics 337 (2001) 23±38

www.elsevier.com/locate/tecto

Dynamic restoration of pro®les across diapiric salt structures: numerical approach and its applications Alik T. Ismail-Zadeh a,b,1,*, Christopher J. Talbot b, Yuri A. Volozh c a

International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Warshavskoye shosse 79, kor. 2, Moscow 113556, Russian Federation b Institute of Earth Sciences, Uppsala University, VillavaÈgen 16, SE-752 36 Uppsala, Sweden c Institute of Geology, Russian Academy of Sciences, Pyzhevsky per. 7, Moscow 109017, Russian Federation Received 5 December 2000; accepted 25 April 2001

Abstract The backstripping method that is widely used in basin analysis sometimes fails for salt-bearing basins because the highly mobile and buoyant salt deforms its sedimentary overburden. We present a numerical approach for 2-D dynamic restoration of cross-sections through successive earlier depositional stages. The approach is based on a solution of the inverse problem of the gravitational (Rayleigh±Taylor) instability and combines the Galerkin-spline ®nite-element method with interface tracking and a backstripping method. Our model interprets basin pro®les as multiple layers of viscous ¯uids with various densities and viscosities. The evolution of salt structures is modelled backward in time by removing successively younger layers and restoring older layers and any diapirs to the stage they were likely to have been. We test the sensitivity of the restoration technique to small variations in density of the layers at different stages in the evolution of diapiric structures. The applicability of the technique was demonstrated by reconstructions of upbuilt and downbuilt diapirs. The technique is used to restore a depthconverted seismic cross-section through the south-eastern part of the Pricaspian salt basin. Mature salt diapirs in the section are shown to have been downbuilt from a salt layer with an initially uniform thickness as a result of differential sedimentary loading until the end of the Triassic before one of the diapirs was buried and actively upbuilt. The numerical approach is well suited for restoration of cross-sections with ductile overburdens, but despite limitations can be developed to 3-D restorations and other rheologies. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Rayleigh±Taylor instability; 2-D backward modelling; Pricaspian basin; salt piercement

1. Introduction

* Corresponding author. Address: International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Warshavskoye shosse 79, kor. 2, Moscow 113556, Russian Federation. E-mail addresses: [email protected] (A.T. Ismail-Zadeh), [email protected] (A.T. Ismail-Zadeh). 1 Present address: Geophysikalisches Institut, UniversitaÈt Karlsruhe, Hertzstr. 16, Geb. 6.42, Karlsruhe 76187, Germany.

Salt tectonics has attracted considerable attention in the last few decades because of its importance in controlling the location of hydrocarbon reserves. Salt diapirs possess physical properties that lead to natural traps of hydrocarbons and allow engineering of cavities for storage of fuels or waste. The evolution of salt structures have been studied numerically by many researchers for the last two

0040-1951/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0040-195 1(01)00111-1

24

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

decades (e.g. Woidt, 1978; Schmeling, 1987; RoÈmer and Neugebauer, 1991; Zaleski and Julien, 1992; Poliakov et al., 1993; Podladchikov et al., 1993; van Keken et al., 1993; Daudre and Cloetingh, 1994; Mello and Henderson, 1997; Naimark et al., 1998). Among them, these models analysed how salt structures related to thickness and viscosity contrasts between the salt and overburden rocks, the effects of local erosion and sedimentation, extension and compression, differential loading and different rheological properties. All these models involved the forward evolution of salt structures toward increasing maturity. To understand the history of deposition, erosion and deformation in sedimentary basins, tools are needed to restore or reconstruct the basin evolution. The backstripping method (Steckler and Watts, 1978) serves as a powerful tool to reconstruct sedimentary layers to their earlier stages. This method is based on stripping off the upper layers step by step and decompacting the lower sediments. Cross-section restoration techniques were initially developed and applied to compressional regions (Dahlstrom, 1968) and later to extensional provinces (Gibbs, 1983). Section restoration enables geoscientists to choose between alternative interpretations of subsurface geometries, analyse the burial history, and determine the evolution of sedimentary basins. At the same time, the methods are limited in the study of salt basins, where it is dif®cult to distinguish the motion of salt and deformation in the overburden during deposition of sediments. A numerical technique for restoration of saltrelated structures was presented by Schultz-Ela (1992) to constrain relevant deformation in extensional terranes. Restoration of a cross-section to the time of deposition of the downward older horizon involved the following successive operations: stripping the upper layer, unfolding deeper deformed horizons, rotation and translation of faulted sedimentary blocks, and decompaction of sediments. Rowan (1993) proposed a modi®cation that included isostatic compensation and/or thermal subsidence. However, these restoration models did not consider motion in the salt itself. In this paper, we present a numerical approach to 2D dynamic section restoration. The approach is based on solving an inverse problem of the gravitational

instability (Ismail-Zadeh et al., 2000a). We employ (1) the Eulerian ®nite element method with interface tracking (Naimark et al., 1998) to determine the deformation of salt and its overburden and (2) the backstripping method (Steckler and Watts, 1978) to strip successive overburden layers, decompact older layers, and `unstrain' the salt structures. We analyse a few synthetic examples to show how the restorations depend on slight changes in density of layers in the sedimentary overburden and on the maturity of salt structures. We demonstrate our approach by restoring two types of salt structures: upbuilt and downbuilt diapirs. The ®rst model represents active diapir growth (i.e. upbuilding), during which all the overburden is prekinematic and was in place before halokinesis began. The second model represents passive diapir growth (i.e. downbuilding), in which accumulation of a thickening overburden depresses the underlying salt while the top of salt in the diapir stays near the depositional (top) surface. The restoration technique is applied to restore a depth-converted seismic cross-section through the south-eastern part of the Pricaspian salt basin on the basis of current hypotheses of the regional geology and tectonics. For the restoration of geological cross-sections, we make the following simplifying assumptions. The restoration is restricted to 2-D deformation of the sections, and there is no out of section deformation. A Newtonian rheology was assumed both for the salt and its overburden. We recognise that most salt overburdens display more complex rheology (e.g. Weijermars et al., 1993). However, there is also evidence that overburden can be considered highly viscous over long time periods. Natural overburdens in a few salt basins exhibit ductile behaviour, for example, Central Iran (Jackson et al., 1995) and the Afgan±Tagik salt basin. Sometimes, overburdens consist of shale and/or impure salt, as in the Pricaspian basin, where pure Kungurian (260 Ma) salt is overlain by impure Kazanian (258 Ma) salt (Volozh et al., 1997, 2001). We do not consider lateral tectonic and thermal effects, although both may play a signi®cant role in halokinesis (Jackson and Talbot, 1994). No fault can be taken into consideration in the restorations, because our approach is based on mechanics of continua rather than fracture mechanics.

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

25

described by the Stokes equation (Landau and Lifshitz, 1987): 27P 1 7 X t 1 rge ˆ 0;

…1†

where P is the pressure, t the deviatoric stress tensor, r the density, and g the acceleration due to gravity, and e the unit vector (0,1). The components of the stress tensor are represented as: ! 2vj 2vi ; i; j ˆ 1; 2; …2† tij ˆ m 1 2xj 2xi

Fig. 1. Geometry of the model region with boundary conditions as indicated. The region is divided into two subregions: V 1 and V 2 by curve L.

2. Model concept We consider a two-dimensional region V …0 # x # L; 2H # z # 0† ®lled by viscous ¯uid and divided by curves La, a ˆ 1; 2; ¼; A (horizons in the geological sense) into several subregions in order to model a layered system of sedimentary cover (Fig. 1). For simplicity in what follows, we consider only one curve, L, and two subregions, V 1 and V 2, though the number of curves can be arbitrary. Physical properties (density and viscosity) can vary from layer to layer in the model. We represent density r (x,z,t) and viscosity m (x,z,t) as sums of two functions:

r…x; z; t† ˆ r 1 …x; z; t† 1 r2 …x; z; t†; m…x; z; t† ˆ m1 …x; z; t† 1 m2 …x; z; t†; where r 1(x,z,t) and m 1(x,z,t) have the ®rst and second continuous derivatives, whereas r 2(x,z,t) and m 2(x,z,t) are represented in the following form: ( 1 ( 1 r2 in V 1 ; m2 in V 1 ; ; r2 ˆ m2 ˆ 2 r2 in V 2 ; m22 in V 2 ; where r12 ; r22 ; m12 ; and m22 are constants. The dynamics of the salt and its overburden are

where n ˆ …v1 ; v2 † is the velocity, x1 ˆ x; x2 ˆ z: The ¯ow is assumed to be incompressible, and hence, we introduce stream function c and represent velocity as v 1 ˆ …2c=2z†; v 2 ˆ 2…2c=2x†: Differentiating the x component of Eq. (1) with respect to z, the z component of the same equation with respect to x, and subtracting one result from the other, we obtain the equation of motion for our layered model: ! ! 22 22 c 22 22 22 c 22 c 4 1 m 2 2 m 2 2x 2z 2x 2z 2z2 2x 2z 2 2x 2 ˆ 2g

2r : 2x

(3)

The advection of density and viscosity is described by: 2r 2c 2r 2c 2r ˆ 2 ; 2t 2x 2z 2z 2x 2m 2c 2m 2c 2m ˆ 2 : 2t 2x 2z 2z 2x

…4†

Curve L: x ˆ X…t†; z ˆ Z…t† satis®es the differential equations: dX 2c dZ 2c ˆ ˆ2 ; : dt dt 2z 2x

…5†

The horizontal boundaries of the model region satisfy the impenetrability condition …v 2 ˆ 0†: Freeslip …t12 ˆ 0† and no-slip …v1 ˆ 0† conditions are assumed for the upper and lower boundaries, respectively. We use symmetry conditions for the vertical boundaries (Fig. 1). In addition, we assume that 2r=2x ˆ 2m=2x ˆ 0 at x ˆ 0 and x ˆ L: This condition is

26

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

necessary to meet the requirements for symmetry and smoothness of r and m in x. The initial conditions take the form r…x; z; t ˆ tin † ˆ r0 …x; z†; m…x; z; t ˆ tin † ˆ m0 …x; z†; X…q† ˆ X 0 …q†; and Z…q† ˆ Z 0 …q†; where r 0(x,z) and m 0(x,z) are the initial distribution of density and viscosity, and X 0 and Z 0 are points on the curve L 0 at t ˆ tin : Thus, the problem is reduced to the computation of the stream function c (x,z,t), density r (x,z,t), viscosity m (x,z,t) and the family of curves satisfying Eqs. (3)± (5), the boundary and initial conditions. 3. Restoration technique In this section, we describe the numerical technique for solving the problems of the gravitational instability backward in time. The basic approach is as follows. Small perturbations of interfaces at which density inversions occur in a layered structure lead to gravity driven deformation. Naimark (1987) has showed that the relevant mathematical problem possesses a unique solution at time interval t in # t # t1 : Hence, functions r0 …x; z† and m0 …x; z†, de®ning density and viscosity of the structure at t ˆ tin ; are transformed into functions r…x; z; t ˆ t1 † and m…x; z; t ˆ t1 †; respectively. In the absence of the dissipation term in the governing equations of motion (3)± (5), the change of the functions (density and viscosity) is a one-to-one continuous transformation, i.e. the functions can be reconstructed backward in time to the initial functions r 0 and m 0. The numerical solution of the problem is based on the Galerkin method applied to Eq. (3) and on the fourth-order Runge±Kutta method applied to Eqs. (4) and (5). The method is described in detail by Naimark et al. (1998). By using the numerical method, we construct a full solution to the timedependent problem for the initial and boundary conditions and given physical parameters. We refer to the model in a Eulerian reference frame; as a result, large displacements are examined without rearrangements of approximating elements. We divide region V into rectangular elements: N 2 1 in x and M 2 1 in z that corresponds to M £ N grid. Approximations to unknown functions c , r 1, and m 1 are represented as linear combinations of basis bicubic splines with unknown coef®cients (here and

below, we assume summation over repeated subscripts taking on the following values: i; k; m ˆ 1; ¼; M; j; l; n ˆ 1; ¼N) b ij …x; z†; c ˆ c ij …t†spij …x; z†; r1 ˆ rij …t†sp b ij …x; z†; m1 ˆ mij …t†sp b ij …x; z† are basis bicubic splines where spij …x; z† and sp satisfying the required boundary conditions. These splines are constructed from basis cubic splines si b ij …x; z† ˆ s^i …x†^sj …z†: and s^i : spij …x; z† ˆ si …x†sj …z† and sp Basis cubic and bicubic splines used here are described by Ismail-Zadeh et al. (1996). Curve L is approximated by a polygon whose vertices have coordinates Xb …t†; Zb …t†; b ˆ 1; ¼B. Applying the Galerkin method and substituting the spline representations into the resulting equations, we derive sets of linear algebraic equations for unknowns c ij and ordinary differential equations for r ij ; mij ; X ˆ Xb …t†; and Z ˆ Zb …t†:

c ij …t†Cijkl ˆ rij …t†Fijkl 1 C kl …t†; drij G ˆ rij …t†Eijkl ; dt ijkl dsj …z† dX ; ˆ c ij …t†si …x† dz dt

…6†

dmij G ˆ mij …t†Eijkl ; …7† dt ijkl dZ ds …x† ˆ 2c ij …t† i sj …z†: dt dx …8†

The forms Cijkl, Fijkl, C kl, Gijkl, and Eijkl represent integrals of products of splines and their derivatives over the model region and are given in Appendix A. To solve the system of equations forward in time, we have to ®nd the unknowns r ij(t), m ij(t), c ij(t), Xb (t), and Z b …t† at t ˆ t~ assuming that the unknowns have been computed at t ˆ t, where Dt ˆ t~ 2 t . 0. Here, we describe the procedure of the numerical solution of Eqs. (6)±(8) backward in time, i.e. for Dt , 0. The time interval is divided by points ts, s ˆ 0; 1; 2; ¼S into S subintervals, and t0 ˆ tin . t1 . ¼ . tS21 . tS ˆ 0. Initial values r ij(t0) and mij …t0 † are derived from the conditions r1 …x; z; tin † ˆ rij …tin †^si …x†^sj …z†; m1 …x; z; tin † ˆ mij …tin †^si …x†^sj …z†: Assuming that the unknowns have been computed at t ˆ ts , t0 ; we can ®nd the matrix Cijkl and the right-hand side of Eq. (6) and solve this set of linear equations for c ij. Using the values so found, we can compute the righthand sides of sets of ordinary differential equations

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

27

Fig. 2. Positions of the interface between two layers in forward (FWM) and backward modelling (BWM). Initial position in FWM (solid line), ®nal position in FWM at t ˆ 1464 (bold), restored interface in BWM for three values of density contrast, Dr : 0.3 (dashed); 0.31 (dot-dashed), and 0.29 (dotted).

in Eqs. (7) and (8). We solve these sets by the fourthorder Runge±Kutta method considering Dt s ˆ ts 2 ts11 to be negative. In this way, we ®nd r ij, m ij, Xb , and Zb at the next time step. Time step Dts is chosen automatically to satisfy the condition that the maximum displacement of material points must not exceed a given value Dl that is suf®ciently small to minimise errors. To choose Dts ; we ®nd vmax ˆ max‰…2c=2z† 2 1 …2c=2x† 2 Š1=2 over grid nodes and specify abs…Dts † ˆ Dl=vmax : At each time interval, the restoration procedure consists of two principal steps: (i) computation of deformations of the layers (solving Eqs. (6)±(8) backward in time) and (ii) application of backstripping technique to the layered system. The latter step means the following. A thin sublayer of uniform thickness is stripped off the top layer of sediments. The thickness of this sublayer, h…ˆ va tc † to be stripped is found from the average rate of deposition, va …ˆ hmax ts21 †, where tc is a time step, and ts and

hmax are the age and maximum thickness of the top layer, respectively. On the removal of the upper layer, underlying sequences are decompacted: the amount of vertical expansion of any sequence at a given location is a function of its lithology and the change in overburden thickness (Sclater and Christie, 1980). Layers of overburden are compactible, whereas salt is assumed to be noncompactible. The new top horizon is then restored to a horizontal line (the assumed depositional con®guration) in a series of time steps. 4. Veri®cation of the technique The restoration procedure was tested by two synthetic examples. First, we developed a model (Fig. 2) consisting of two layers with dimensionless density contrast Dr ˆ r1 2 r2 ˆ 0:3 and viscosity ratio n ˆ m1 =m2 ˆ 10; subscripts 1 and 2 refer to the

28

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

upper and lower layers respectively. The aspect ratio ( ˆ model width/model height) was assumed to be 3. The interface between layers (solid line in Fig. 2) was disturbed initially by a two-peak perturbation. The evolution of the layered structure was modeled forward in time. At dimensionless time t ˆ 1464 (500 computational time steps, CTS) the calculations were interrupted. The bold line in Fig. 2 (upper panel) presents the ®nal position of the interface in forward modelling. To restore the interface to its initial position (at t ˆ 0), we used the restoration technique (with no backstripping) in the cases when the density contrast Dr was equal to (i) 0.3, (ii) 0.29, and (iii) 0.31. The ®nal positions of the restored curves for the three cases are presented in Fig. 2. These curves are almost invisible in the upper panel due to their proximity to the curve of the initial perturbation. The bottom panel of Fig. 2 presents an enlarged view of the curves. In the ®rst case, the restored curve (dashed line) matches the initial interface almost exactly. A very small discrepancies (less than 0.1%) between the two curves is the result of computational errors accumulated during 1000 CTS (500 steps in forward modelling 1 500 steps in backward modelling). In cases (ii) and (iii), slight changes in density contrast lead to slight de¯ections of restored interfaces from the position of the initial interface. The test shows a rather high accuracy of the restoration procedure. The second test was an extension of forward computations in the previous model in order to simulate an evolution of a mature diapir. Fig. 3 presents diapir evolution (solid line) for six successive time steps. The computations were stopped at t ˆ 3712 (2500 CTS). We can see the complex shape of the resultant diapir that evolve from an initial perturbation with two peaks. The restoration technique was applied to the ®nal position of the interface between the two layers. The discrepancies between the interface (solid lines in Fig. 3) and its restoration (dashed lines in Fig. 3) at respective times become suf®ciently visible as the result of error accumulation during the computations (after about 4000 CTS ˆ 2500 CTS in forward 1 1500 CTS in backward modelling). If the model is computed for a very long time, the light layer will spread into a horizontal layer near the surface, and practical restoration becomes impossible.

5. Restoration of demonstration models We restore the evolution of two cases of halokinesis. Model adg represents an active piercement of salt, during which a postdepositional diapir grows through prekinematic overburden (Jackson and Talbot, 1994). As the diapir rises, its base remains at a constant depth below the sedimentary surface, while its crest rises toward the surface. Model pdg represents a passive piercement of salt, which describes syndepositional diapir growth: the diapir increases in relief by growing downward relative to the ®xed depositional surface. The base of the diapir subsides, along with surrounding strata, as the salt-withdrawal basin ®lls with sediments. The diapir's crest remains close to the depositional surface although a thin temporary roof may be thickened by sedimentation or thinned by erosion. In both models, salt is assumed to be incompressible and to have a density of 2:2 £ 103 kg m 23. The density of overburden r ov(z) increases exponentially with depth z (Biot and OdeÂ, 1965):

rov …z† ˆ ra 2 rb exp…2az† 2 0:5z exp…2gz†; where r a ˆ 2.5 £ 10 3 kg m 23, r b ˆ 0.29 £ 103 kg m 23, a ˆ 5:9 £ 1024 m21 ; and g ˆ 1:6 £ 1023 m21 . The density of each layer is found by averaging its `density' over grid points. During the computation of deformation of the layered structure, the densities of the layers are not changed. We assume viscosities of 10 20 Pa s for the overburden and 10 18 Pa s for salt. The square model box is 5 km long and 5 km deep. This box is divided into 49 £ 47 rectangular elements in the x and z directions, respectively. The time scale used is tp ˆ mp =…rp gH† ˆ 266:7 yr; where mp ˆ 1018 Pa s; rp ˆ 2:5 £ 103 kg m23 ; and H ˆ 5 km: 5.1. Model adg Initially, we developed a forward model of upbuilt diapir to verify the validity of our restoration procedure. In the forward model, a basal salt layer, 1 km thick, was overlain by a prekinematic overburden 4 km thick consisting of four layers, having different physical properties. The interface is initially perturbed by a peak of cosine shape. The structure is gravitationally unstable, and a small perturbation of the salt/ sediment interface results in diapir growth. The ®nal position of the interfaces between sedimentary layers

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

29

Fig. 3. Shapes of diapir in forward modelling (solid lines) and of restored diapir in backward modelling (dashed lines) at successive times indicated (tf =tb is dimensionless time in forward/backward modelling).

30

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Fig. 4. Dynamic restorations of an upbuilt diapir to times indicated (model adg). In each panel, salt is white, and layers of ductile overburden are marked by shades of grey. Vertical and horizontal scales are dimensionless.

in the forward model was used as an initial position of the interfaces for a backward model. Fig. 4 shows the results of restoration of the upbuilt diapir (lower white layer) from the position at t ˆ 0 to its position for time t ˆ 8:3 Ma: 5.2. Model pdg As in the case of the previous model, to verify our restoration procedure, we initially developed a forward model of downbuilt diapir. We considered a three-layered model consisting of salt (lowest layer)

downbuilt by sediments (middle layer), and a layer above a depositional surface (uppermost layer of air or water). The density of the uppermost layer was assumed to be zero, and its viscosity was several orders of magnitude lower than the viscosity of the overburden. Obviously, the air/water layer in the model is only an approximation to nature; however, it has been shown that such approximation works well enough in many applications (Naimark et al., 1998). We assumed that fast erosion and redeposition maintained the surface topography ¯at over the diapir. During the modelling, horizontal sedimentary layers

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

31

Fig. 5. Dynamic restorations of a downbuilt diapir to times indicated (model pdg). Shading and scaling is as in Fig. 4.

of uniform thickness (referred to as non-kinematic layers) were introduced above the sedimentation surface: the sediment deposition was considered to be rapid compared to diapir downbuilding. At the same time, we changed the density and thickness of each sedimentary layer below the non-kinematic layer according to the model of compaction (Biot and OdeÂ, 1965; Sclater and Christie, 1980). The forward modelling was continued until the sedimentary overburden consisted of ®ve layers. The ®nal position of forward modelling of the multi-layered structure was taken as an initial con®guration for dynamic restoration. The model of a salt diapir (white) downbuilt by the ®ve layers of overburden (grey shadings) was

restored from the present time to its position for time t ˆ 23:9 Ma (Fig. 5). Let us describe the restoration procedure in model pdg. The layered structure was numerically restored to a position where an upper layer becomes nearly horizontal as it was during the introduction of the non-kinematic layer in forward modelling. The upper layer was then stripped off, and the underlying layers were decompacted according to a standard porosity- and density-depth exponential law. After stripping successive upper layer from the model, the values of `depth' and density for each grid point were computed within the remaining underlying layers. The procedure was applied successively to the layered system until a single layer of overburden was left.

32

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Fig. 6. Map of the Pricaspian basin showing salt structures, major structural features, and regions I±VI distinguished by morphology of salt structures: (I) equidimensional salt pillows; (II) wall-like salt pillows; (III) salt anticlines and turtle-backs; (IV) salt domes, some have surfaced and extruded or recycled; (V) salt diapirs and walls sometimes complicated by salt overhangs; (VI) salt domes and walls where halokinesis involved both Kungurian and Kazanian salt. (1) pre-Jurassic erosional boundary of salt deposits; (2) pre-Jurassic thrust front; (3) interfaces between regions; (4) grey/white ˆ areas of salt thicker/thinner than the initial salt deposits; (5) overhangs; (6) undisturbed bedded salt on platform; (7) outline of the Caspian Sea.

For each demostration model (adg and pdg), we compared the layered structure predicted by the restoration with the initial structure of the forward model for the same time interval. Very small computational errors (about 0.1%) were observed in these comparisons.

6. Restoration of cross-section in the Pricaspian salt basin 6.1. Interplay between sedimentation and halokinesis The western literature often uses the adjectives `Pre-Caspian' (earlier than Caspian) or `Peri-Caspian'

(peripheral to Caspian) to refer to the basin that is approximately 600 km across from west to east and is underlain by salt at the northern end of the Caspian Sea at the south-eastern portion of the East-European platform (Fig. 6). To avoid such ambiguities in time or space, we use here the Russian `Pricaspian' which implies an area adjoining the Caspian. The sedimentary basin in®ll is divided into three major sequences: subsalt, salt, and salt overburden (Konishev and Volozh, 1990; Brunet et al., 1999). The subsalt sequence represents a number of stratigraphic horizons separated by unconformities, namely the Riphean±Lower Palaeozoic, Devonian±Lower Carboniferous, and Middle Carboniferous±Lower

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Permian horizons. The subsalt sequence has a fairly complex structure, although it is dominated overall by two lithologies (carbonate and terrigenous) which pass into each other both laterally and vertically. The salt sequence in the basin is predominantly of Kungurian time (260±258 Ma). The Kungurian salt was overlain by Kazanian (258±252 Ma) salt in the centre of the basin, whereas it was downbuilt along the eastern and south-eastern margins by Kazanian clastic sediments prograding from the rising Urals (Volozh et al., 2001). The Permian salt developed a variety of salt structures, such as salt pillows, 2-D diapiric walls and 3-D stocks complicated by extensive overhangs and giant salt massifs (Volozh et al., 1994, 1997). The salt overburden consists predominantly of terrigenous Upper Permian, Triassic, Jurassic, Cretaceous, Palaeogene, and Neogene sediments. Carbonate formations occur in the Upper Cretaceous, Upper Jurassic, and, in the western part of the basin, Middle Triassic. The overburden is divided into three structural levels by gentle unconformities: Upper Permian±Triassic, Jurassic±Miocene, and Pliocene± Quaternary. Structures in the overburden are determined by the location of salt structures. Salt-withdrawal basins are clearly visible in the Upper Permian± Triassic sediments. Jurassic±Palaeogene sediments lie unconformably on the Upper Permian±Triassic sequence and are less intensively deformed than the underlying unit. Pliocene sediments in®ll prePliocene erosional features. Permian salt was extruded as allochthonous sheets and/or was redeposited as salt beds from Triassic to Jurassic and, probably, Cretaceous times (Volozh et al., 1989, 2001). A number of regions of the Pricaspian salt basin can be distinguished on the basis of the morphology of their salt structures (Volozh et al., 1989; Konishev and Volozh, 1989; Volozh et al., 2001). Salt structures are immature in a north-eastern region (I in Fig. 6) (over a basement ridge). Here salt pillows are buried conformably by Upper Permian±Triassic deposits. Region II along the north-western ¯ank of the basin is characterised by several parallel symmetrical salt walls extending over hundreds of kilometre long. The height of the salt walls ranges from 3 to 5 km and increases progressively towards the centre of the basin. Region III along the south-east margin of the Pricaspian basin is marked by series of linear elongate

33

asymmetric salt walls with a gently dipping outer ¯anks and steeply dipping basinward ¯anks. The salt domes in region IV are mature and often complicated by overhangs and sometimes by salt-tongues or even canopies. Dimension of the overhangs vary widely (length from 1.5±2 to 25±30 km, width from 1.5±2 to 4±6 km). A region of mature salt structures (V) is situated to the west and north-east of region IV. A region of large immature salt structures (VI) occupies the central portion of the basin. A signi®cant feature of the Pricaspian basin is that the majority of salt structures with a wide variety of shapes and amplitudes are joined into chains (Volozh et al., 2001). 6.2. Restoration model To apply our restoration technique to geological structures, we chose a depth-converted part of a NW±SE seismic pro®le through the south-eastern border of the Pricaspian basin (region III in Fig. 6) where salt walls and their intervening salt-withdrawal basin have 2-D geometries (Fig. 7). The salt-withdrawal basin between two diapiric salt walls was in®lled by Upper Permian clastic sediments which are unconformably overlain by thick Triassic±Cretaceous and thin Paleogene±Quaternary terrigeneous-carbonate deposits. The maximum thickness of the overburden is about 4.5 km. We introduce 21 layers in the restoration model: a basal salt layer, 19 layers of sedimentary overburden, and a layer above the current depositional surface. A few faults in the overburden were removed from the geological section (Fig. 7), because (i) our viscous model cannot handle them, and (ii) these faults have a little or no effects on restoration of the layers. We used the relationship between the depth, h (km), and density of the overburden, r ov (g cm 23), suggested by Nevolin et al. (1977) for the Pricaspian basin: 8 > 1:75 1 0:45h 2 0:01h2 ; if h , 1:1; > < rov ˆ 2:02 1 0:21h 2 0:0175h2 ; if 1:1 # h , 6; > > : 2:62 1 0:05h; if 6 # h , 10: The average density and viscosity of salt were assumed to be 2:24 £ 103 kg m 23 and 10 17 Pa s, respectively (Nevolin et al., 1977). Terrigeneous clastic and carbonate sediments such as those in the Pricaspian region typically have a higher effective

34

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Fig. 7. (A) Seismic pro®le through the south-eastern part of the Pricaspian basin (see Fig. 6. for location) and (B) its depth-conversion with salt related structures as indicated: (1) salt dome; (2) residual salt high; (3±8) sedimentary layers of salt overburden ((3) Kazanian, Upper Permian (258±252 Ma); (4) Tatarian, Upper Permian (252±248 Ma); (5) Upper Triassic (246±242 Ma); (6) Middle to Lower Triassic (242±205); (7) Middle Jurassic to Lower Cretaceous (170±98 Ma); and (8) Upper Cretaceous to Quaternary (98±0 Ma)).

viscosity than evaporites. We treated the viscosity ratio between salt and overburden as the leastknown parameter. As a modelling time step, t, is proportional to viscosity value, m (i.e. t ˆ m=‰rp gHŠ), the viscosity of the overburden was used as a tuning

parameter of our model in order to keep time of model evolution close to observed geological times. The overburden was taken as 10 to 10 3 more viscous than salt. Fig. 8 shows the reconstructed history of the

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38 Fig. 8. Restoration model of the depth-converted seismic pro®le from the south-eastern part of the Pricaspian basin shown in Fig. 7B. Triangle indicates a depocentre's position.

35

36

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

cross-section predicted by the suggested method. Salt was clearly mobile along the whole length of the section, as shown by the lateral variation in salt and clastic bed thicknesses. Despite any possible appearance to the contrary, the area of salt in the restored pro®les remains constant. We applied a minimalist approach to unconformities (see panels at 170, 215, and 228 Ma in Fig. 8) to avoid uncertainties concerning any volume of sediments eroded. Hence, we ignore the possibility of erosion by removing only the overburden actually present. Although the restoration was carried out backward in time, the results can be discussed in terms of forward evolution. Differential loading by sediments supplied from Urals initiated salt ¯ow in Early Kazanian time. Once the salt movement began, the dominant process controlling the dynamics of the salt was local downbuilding of its overburden. It is clearly seen that a depocentre moved to the southeast till the Middle Triassic (,238 Ma). Salt was squeezed from beneath a deepening salt-withdrawal basin and driven into the diapirs as the remaining portions of the mother salt layer developed some relief by differential thinning and reorganisation. A turtleback structure was developed at depth as the salt withdrew from the source layer. The crest of the diapirs remained near the surface till end-Triassic. It is likely that some salt was extruded and eroded or dissolved. However, we avoided uncertainties concerning volume of the salt eroded or dissolved by keeping the area of salt constant. After the hiatus (,205± 170 Ma) a distinguishing feature of salt dynamics became an active piercement of the diapir at the right-hand end of the section (in the south-east) through the prekinematic overburden accumulated from the Upper Jurassic to Paleogene. The diapir at the left-hand end of the section (in the north-west) stagnated as it was isolated by a primary weld. The evolution of the salt and its overburden recovered is in agreement with geological and tectonic hypotheses on the development of the region. 7. Conclusions We have presented a numerical method for the dynamic restoration of salt and its overburden back to the salt depositional phase. The technique combines

the Galerkin-spline ®nite-element method with procedures for interface tracking and backstripping. Testing of the restoration technique showed its stability to small changes of initial conditions. We demonstrated the practicality of the suggested technique by reconstructions of the histories of model upbuilt and downbuilt diapirs. The technique was used to restore a depth-converted pro®le through the south-eastern margin of the Pricaspian basin. That pro®le is characterised by mature salt diapirs which were shown to have been downbuilt from an initially continuous salt layer modi®ed by differential sedimentary loading and upbuilt after being buried at ,260 Ma. The power of the method suggested should not be overestimated for it still has many geological limitations. It can only restore 2-D pro®les of salt diapirs with ductile overburdens and cannot handle signi®cant faults. At the same time, the restoration technique can be extended to 3-D problems and more complex rheologies of the sedimentary overburden. Recently, Ismail-Zadeh et al. (2000b), Kaus and Podladchikov (2001) and Korotkii et al. (2001) developed numerical techniques to compute 3-D forward and backward models of Rayleigh±Taylor instability. These techniques can be employed for restorations of three-dimensional diapiric structures. Besides the method's basic limitations, there are also several sources of errors. Pro®le restoration depends strongly on the initial data set, namely, seismic velocities used for depth conversion of seismic pro®les, lithologies speci®ed for decompaction, and ages assigned to the interpreted horizons. The method also depends on the density±velocity and effective viscosity with depth relationships for the speci®c region and site speci®c sediment compaction curves. As presented, to avoid non-unique solutions we ignored loss of overburden by erosion and salt by extrusion and dissolution, although this could be easily added to the models. Despite its limitations this technique provides a quantitative analytical tool to study the interplay between sedimentation and salt mobility. In a more comprehensive sense, dynamic restoration of crosssections provides spatial±temporal information on the formation and evolution of structural traps and a general framework for evaluating hydrocarbon migration pathways.

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Acknowledgements We are very grateful to B. Naimark, H. Schmeling, B. Vendeville, and R. Weijermars for critical comments which signi®cantly improved an early version of the manuscript. We are thankful to A. Poliakov and R. Weinberg for their constructive reviews of the manuscript. The publication was supported by the Royal Swedish Academy of Sciences (KVA 1325) and the International Science and Technology Center (ISTC 1293). A.T.I.-Z. was also supported by the Fellowship of the President of the Russian Federation to Young DSc (99-15-96085), RFBR (99-05-65050) and US National Science Foundation (EAR 9804859). Appendix A. Coef®cients in the discrete governing equations Here, we present coef®cients entering into Eqs. (6) and (7) (see Naimark et al., 1998 for more detail). Coef®cients Cijkl on the left-hand side of Eq. (6) are represented as sums of three terms: Cijkl ˆ 1 2 3 Cijkl 1 Cijkl 1 Cijkl where 1 110 000 220 200 020 Cijkl ˆ mmn …4A110 ikm Bjln 1 Aikm Bjln 2 Aikm Bjln

Apgr ikm ˆ Bpgr jln

Z1 0

37

function (´´´) and the zero-order derivative is the func1 depend on the tion itself. We see that elements Cijkl continuous term m 1, but are independent of the curve 2 3 L. On the other hand, elements Cijkl and Cijkl depend 1 on the curve L and on the constants m2 and m22 , but are independent of the continuous term m 1. Coef®cients Fijkl and C kl on the right-hand side of Eq. (6) are represented in the following form: Z1 Z1 F ijkl ˆ …^si …x††…1† sk …x† dx s^j …z†sl …z† dz; 0

C kl ˆ g…r22 2 r12 †

0

Z L

sk …x†sl …z† dz:

Coef®cients Gijkl and Eijkl entering Eq. (7) are calculated by integrating basic splines and their derivatives: Z1 Z1 Gijkl ˆ s^i …x†^sk …x† dx s^j …z†^sl …z† dz; 0

0

100 100 001 Eijkl ˆ c mn …A001 ikm Bjln 2 Aikm Bjln †; pqr pqr pqr pqr where Aikm and Bjln are obtained from Aikm and Bjln with si …x†, sk …x†, s^m …x†; sj …z†, sl …z†, sl(z), and s^n …z† replaced by s^i …x†; s^k …x†; sm(x), s^j …z†; s^l …z† and sn(z), respectively.

200 220 000 2 A020 ikm Bjln 1 Aikm Bjln †;

References

Z1

Biot, M.A., OdeÂ, H., 1965. Theory of gravity instability with variable overburden and compaction. Geophysics 30, 213±227. Brunet, M.-F., Volozh, Yu.A., Antipov, M.P., Lobkovsky, L.I., 1999. The geodynamic evolution of the Precaspian Basin (Kazakhstan) along a north±south section. Tectonophysics 313, 85±106. Dahlstrom, C.D.A., 1968. Balanced cross sections. Can. J. Earth Sci. 6, 743±757. DaudreÂ, B., Cloetingh, S., 1994. Numerical modelling of salt diapirism: in¯uence of the tectonic regime. Tectonophysics 240, 59± 79. Gibbs, A.D., 1983. Balanced cross-section construction from seismic sections in areas of extensional tectonics. J. Struct. Geol. 5, 153±160. Ismail-Zadeh, A.T., Naimark, B.M., Lobkovsky, L.I., 1996. Hydrodynamic model of sedimentary basin formation based on development and subsequent phase transformation of a magmatic lens in the upper mantle. In: Chowdhury, D.K. (Ed.), Computational Seismology and Geodynamics, vol. 3. American Geophysical Union, Washington, DC, pp. 42±53. Ismail-Zadeh, A.T., Talbot, C.J., Naimark, B.M., 2000a. Reconstructing the evolution of layered geostructures: inverse problem of gravitational instability. In: Keilis-Borok, V.I.,

0

si …x†…p† sk …x†…q† s^m …x†…r† dx;

sj …z†…p† sl …z†…q† s^n …z†…r† dz;

2 Cijkl ˆ m12

3 ˆ m22 Cijkl

Z

Z V1

Z

Z V2

…L…si …x†sj …z†; sk …x†sl …z†† dx dz; …L…si …x†sj …z†; sk …x†s1 …z†† dx dz;

2 2c 2 2w 1 L…c; w† ˆ 4 2x 2z 2x 2z 

22 c 22 c 2 2z 2 2x 2

!

! 22 w 22 w 2 2 ; 2z2 2x

where (´ ´ ´) (p) denotes the derivative of order p of a

38

A.T. Ismail-Zadeh et al. / Tectonophysics 337 (2001) 23±38

Molchan, G.M. (Eds.), Problems in Geodynamics and Seismicity of the Earth. GEOS, Moscow, pp. 52±61. Ismail-Zadeh, A.T., Tsepelev, I.A., Talbot, C., Oster, P., 2000b. A numerical method and parallel algorithm for three-dimensional modeling of salt diapirism. In: Keilis-Borok, V.I., Molchan, G.M. (Eds.), Problems in Dynamics and Seismicity of the Earth. GEOS, Moscow, pp. 62±76. Jackson, M.P.A., Talbot, C.J., 1994. Advances in salt tectonics. In: Hancock, P.L. (Ed.), Continental Deformation. Pergamon Press, Oxford, pp. 159±179. Jackson, M.P.A., Roberts, D.G., Snelson, S., 1995. Salt tectonics Ð a global perspective. AAPG Mem. 65, Tulsa, 454 pp. Kaus, B.J., Podladchikov, Y.Y., 2001. Forward and reverse modeling of the three-dimensional viscous Rayleigh±Taylor instability. Geophys. Res. Lett. 28 (6), 1095±1098. van Keken, P.E., Spiers, C.J., van den Berg, A.P., Muyzert, E.J., 1993. The effective viscosity of rocksalt: implementation of steady-state creep laws in numerical models of salt diapirism. Tectonophysics 225, 457±476. Konishev, V.S., Volozh, Yu.A., 1989. Tectonic zonation of salt structures in areas of halogenesis. Dokl./Proc. Belyaruss. Acad. Sci. 33 (9), 832±835 in Russian. Konishev, V.S., Volozh, Yu.A., 1990. Formation of salt structures at the continental margins. Dokl./Proc. Russ. Acad. Sci. 310 (4), 939±941 in Russian. Korotkii, A.I., Tsepelev, I.A., Ismail-Zadeh, A.T., and Naimark, B.M., 2001. Three-dimensional backward modeling in problems of Rayleigh-Taylor instability. Proc. Ural State Univ. 22, 1-9 in Russian. Landau, L.D., Lifshitz, E.M., 1987. Fluid Mechanics. 2nd ed. Pergamon Press, Oxford. Mello, U.T., Henderson, M.E., 1997. Techniques for including large deformations associated with salt and fault motion in basin modeling. Mar. Petrol. Geol. 14 (5), 551±564. Naimark, B.M., 1987. Existence and uniqueness of the solution of the Rayleigh±Taylor problem. Computat. Seismol. 18, 32±41. Naimark, B.M., Ismail-Zadeh, A.T., Jacoby, W.R., 1998. Numerical approach to problems of gravitational instability of geostructures with advected material boundaries. Geophys. J. Int. 134, 473±483. Nevolin, N.B., Kunin, N.Ya., Andreev, A.P., Volozh, Yu.A., 1977. Based on geophysical data geological structure and oil±gas content of intracratonic salt basins. Nedra, Moscow in Russian.

Podladchikov, Yu., Talbot, C., Poliakov, A., 1993. Numerical models of complex diapirs. Tectonophysics 228, 189±198. Poliakov, A., van Balen, R., Podladchikov, Yu., Daudre, B., Cloetingh, S., Talbot, C., 1993. Numerical analysis of how sedimentation and redistribution of sur®cial sediments affects salt diapirism. Tectonophysics 226, 199±216. Rowan, M.G., 1993. A systematic technique for the sequential restoration of salt structures. Tectonophysics 228, 331±348. RoÈmer, M.-M., Neugebauer, H.J., 1991. The salt dome problem: a multilayered approach. J. Geophys. Res. 96, 2389±2396. Schmeling, H., 1987. On the relation between initial conditions and late stages of Rayleigh±Taylor instabilities. Tectonophysics 133, 65±80. Schultz-Ela, D., 1992. Restoration of cross-sections to constrain deformation processes of extensional terranes. Mar. Petrol. Geol. 9, 372±388. Sclater, J.G., Christie, P.A.F., 1980. Continental stretching: an explanation of the post-mid-Cretaceous subsidence of the Central North Sea Basin. J. Geophys. Res. 85, 3711±3739. Steckler, M.S., Watts, A.B., 1978. Subsidence of the Atlantic type continental margin off New York. Earth Planet. Sci. Lett. 42, 1±13. Volozh, Yu.A., Votsalevsky, E.S., Zhivoderov, A.B., Nurbaev, B.O., Pilifosov, V.M., 1989. Problems of oil and foulness of sediments above salt in the Pricaspian depression. Izvest., Kazakh Acad. Sci., Geol. Ser. 4, 3±15 in Russian. Volozh, Yu.A., Groshev, V.G., Sinelnikov, A.V., 1994. The overhangs of the southern Precaspian basin (Kazakhstan): proposals for a genetic classi®cation. Bull. Centre Res. Explor. Elf Aquit. 18, 19±32. Volozh, Yu.A., Volchegurskii, L.F., Groshev, V.G., Shishkina, T.Yu., 1997. Types of salt structures in the Peri-Caspian depression. Geotectonics 31, 204±217. Volozh, Yu.A., Talbot, C.J., Ismail-Zadeh, A.T., 2001. Salt structures and hydrocarbons in the Pricaspian basin. AAPG Bull., submitted for publication. Weijermars, R., Jackson, M.P.A., Vendeville, B., 1993. Rheological and tectonic modeling of salt provinces. Tectonophysics 217, 143±174. Woidt, W.-D., 1978. Finite element calculations applied to salt dome analysis. Tectonophysics 50, 369±386. Zaleski, S., Julien, P., 1992. Numerical simulation of Rayleigh± Taylor instability for single and multiple salt diapirs. Tectonophysics 206, 55±69.