- Email: [email protected]
Finite element models of density instabilities by means of bicubic spline interpolation
176
Physics of the Earth and Planetary Interiors, 21(1980)176—180 © Elsevier Scientific Publishing Company, Amsterdam — Printed in The Netherlands
FINITE ELEMENT MODELS OF DENSITY INSTABILITIES BY MEANS OF BICUBIC SPLINE INTERPOLATION W.-D. WOIDT
1
and HJ.
2
1 Inst itut für Geophysik, TU Braunschweig, Brunswick (Federal Republic of Germany) 2 Institu t für Geophysik, Goethe Universitàt, 6000 Frankfurt (Federal Republic of Germany)
(Received February 14, 1979; accepted for publication May 12, 1979)
Woidt, W.-D. and Neugebauer, H.J., 1980. Finite Element models of density instabilities by means of bicubic spline interpolation. Phys. Earth Planet. Inter., 21: 176—180. Density instabilities are often associated with large deformations in rocks, and result in diapiric structures of different shape and dimension. The Finite Element model presented here is able to follow the large deformations in a simple two-layer model with inverted density stratification. The developing interface instability gives interesting insights into the dependence of the rising diapiric structures on the model geometry, and the viscosity contrast between the two fluids. The development of single diapiric structures, as well as families of diapirs, are studied with the Finite Element model.
1. Introduction A model consisting of two superimposed viscous fluids with inverted density stratification, which develops spontaneously into a series of ridges and domes, has often been used to explain the dynamical behaviour of geological formations like salt domes, salt dome families, gneiss domes and, on a larger scale, continental rift systems and plumes. Such an unstable situation between fluids of different density and viscosity is generally called a Rayleigh—Taylor instability. As the viscosities of the fluids are very high if applied to geological timescales, the convective terms in the Navier—Stokes equation can be neglected, and nonlinearity only appears in the boundary conditions at the interface. For geophysical considerations, the temporal and spatial development of the interface between the two different fluids is of specific interest. Its shape can be studied as a function of the viscosity contrast and the boundary conditions at the interface. Figure 1 gives a short outline of the model geometry and the boundary conditions,
Analytical solutions of this two-dimensional model have been intensively studied in connection with the analysis of salt domes (Danes, 1964; Selig, 1965; Ramberg, 1968) and salt dome families (Selig and Wermund, 1966; Hunsche, 1977). These analytical solutions linearize the nonlinear term of the interface equation, and are therefore limited to infinitesimal disturbances. The linear model predicts an exponential growth of every sinusoidal disturbance in the interface. Disturbances of one specific wavenumber are found to show maximum growth rate, and therefore will determine the spacing of the developing interface structures. Because the linear model is limited to infmitesimal disturbances, it does not allow any assumption about the shape of the growing instability. The results of the analytical analysis can be applied to any artifical initial perturbation of the interface by decomposition into Fourier components. Every fundamental mode develops corresponding to its specific growth rate. After a time-increment, which is limited by the applicability of the linear analysis, the superposition of all independently
177
boundary conditions: z=h
1
:
0 or vI
=
:
Z=~(xt)
zs—h2 :
[ngn]~~
0
,
=
0
0. a’
=
=
0
0
~
=
•
0
0
Fig. 1. Two-dimensional model geometry and the boundary conditions.
growing components shows a wave-front propagation triggered by the initial perturbation. Such models have been successfully used to explain the phenomena of salt dome families, Starting from the results of the linear model, Nasir and Dabousi (1978) constructed a general solution by using’an infinite series of eigenfunctions. The infinite series solution is truncated to obtain the constants of integration from the boundary conditions along the interface, including the nonlinear term. Therefore their model makes it possible to study the shape of the growing instability. The results show the transition from the linear to the nonlinear development of the interface instability, but no diapiric stage. Up to now, diapiric structures have been studied mainly by model experiments using oils and glycerine (e.g., Whitehead and Luther, 1975) or with pasty materials in a centrifuge (e.g., Ramberg, 1967).
2. Numerical models In numerical calculations with the Finite Element method, the authors have determined that the main difficulty lies in following the large deformations in the two fluids. For this purpose, a deformable mesh, which follows the movement of the fluids, has only limited flexibility. The advantage of such a deformable mesh, however, is that the boundary conditions at the interface between different materials can be treated as natural, and therefore are easy to fulfil. The Finite Element model used here is based on the stream function approach, in which the condi-
tion of incompressibility is automatically implemented at the cost of fourth-order derivatives in the differential equation (Woidt, 1979). The interface between the two fluids is represented by discrete points, moving through a mesh, which remains stationary in space. Consequently, the abrupt change of the viscosities at the interface can only be described approximately in this model. Rectangular elements with high-order interpolation functions are used, which have to be of C’-continuity at least. Bicubic Hermitian polynomials (Prenter, 1975) as interpolation functions shows a very good order of accuracy, but unfortunately result in a large overall matrix system due to the sixteen unknowns per element. To follow the interface for as many time-steps as are necessary to reach its fully developed diapiric stage, within the capability of present day computers, requires basically small dimensions in the overall matrix. The number of unknowns can be considerably reduced by using bicubic spUme interpolation functions,which have a local support over sixteen elements and which are continuously differentiable twice (Schultz, 1973). The only unknown in this approach is the stream function itself. A typical mesh size of 16 X 25 rectangular elements results in 1768 unknowns and a semiband width of 72, using the bicubic Hermitian polynomials as interpolation functions in the Finite Element calculations. In comparison, the more complicated numerical algorithm, which makes use of bicubic spline interpolation functions, results in an overall matrix with 360 unknowns, and a semiband width of 49 for the same mesh.
178
The time integration scheme by which the interface is moved into its new position is a fourth-order Runge—Kutta or Adam—Bashford process. 40—100 time-steps are necessary to reach the fully developed diapiric stage. Lower-order schemes give wrong results because of the exponential or even over-exponential growth of the disturbance in the beginning of the dome’s development.
I
t~ 22.6
.~ “
20 ,
h1 18
~-.-
3. Results and discussion The first two plots (Figs. 2 and 3) show diapiric structures, which evolve from an initially sinusoidal perturbation of the interface. The only difference in the parameters of both models is in the ratio of viscosities, which is 100 : 1 in Fig. 2 and 1 :100 in Fig. 3. As it can be seen in the plots, the shape of the growing interface instability is very much dependent on the ratio of viscosities. Following the top of the rising diapir in Fig. 3, an over-exponential growth can be derived during the nonlinear stage, which is characteristic for viscosity ratios clearly less than unity. The thickness of the heavier upper layer is seven times that of the underlying lighter fluid. Because of this, the rising instability can develop its different stages without being influenced too early by the upper boundary. Such an early influence forces the upwelling lighter fluid to flow aside, forming a mushroom-shaped structure. The strong influence of this boundary on the development of the t~ ~ h
_________
.~—
8 _________
P
‘~
I Fig. 3. As for Fig. 2, except that the ratio of viscosities ~I/M2 = 0.01, and A = 2.96 h2. ,
instability can be studied in Fig. 4, where the thickness of both layers is about equal, and the ratio ofviscosities is 100 : 1. Compared with Fig. 2, the stage where the rising fluid forms a spherical bulb in Fig. 4 is suppressed by the upper boundary condition. In each case the wavenumber of the initial disturbance shows maximum growth as predicted by the linear analytical solution. Good agreement between the analytical solution and this numerical model has always been found for the linear stage of the interface instability. Figures 5 and 6 show diapiric structures where the initial triangular-shaped disturbance triggers the evolution of new diapirs From model experiments with oils it is known that viscosity ratios near to unity form the most fascinating instabilities (Heye,
~
1 h2:
17
13
~2
Fig 2 Model diapir which develops from an initially sinu soidal perturbation h1/h2 ~ = 100 and 3.A Free-slip = 13 h2 conditions on the boundary fixed (zero velocThe densityhold contrast forupper all models is 0.1and g cm ities) conditions on the lower surface. Time can be derived from the formula t = t*M2/h 2 (P2 in poises, h2 (cm) and t (s)). can be taken directly from the plots (see text).
1978). Therefore, a viscosity ratio of unity is chosen hi h2
__
________________________________________ .
.
—
—
.
—
Fig. 4. As for Fig. 2, except that the ratio of thicknesses h1/h2 = 0.92, and therefore A 5.36 h2.
179
boundary are exchanged, as has been done in Fig. 6, the dominant wavelength becomes much greater, and
~‘
h
2 1.0- 1.2- 1.32
Fig. 5. Model diapirs which are triggered by a triangular initial perturbation (width 6h2). Model parameters are ~ /~2= 1, hi/h2 = 7, and the length of the modelL = 48h2.
for these numerical calculations. The two examples differ only in the conditions on the upper and lower boundaries, which are free-slip on the upper boundary and non-slip on the lower boundary in Fig. 5, and vice versa in Fig. 6. The wavelength with maximum growth rate, as predicted by the linear solution, is expected, to dominate the distance between the different generations of diapirs,which is 3.7 times the thickness of the thin layer for Fig. 5. This theoretical prediction is in good agreement with the numerical calculations for the first generation of diapirs, far first as the early nonlinear stage is concerned. Later,asthe generation begins to grow into the direction of the mother diapir, and the distance between them becomes shorter. Finally, material streams into the second rim syncline, thereby distorting the first generation diapirs into a direction away from the mother diapir. It is interesting to note that the distance between first and second generation diapirs is much greater than that theoretically expected. This can be understood if one considers that the lighter fluid layer is too thin to feed a new diapiric structure every 3.7 h2, and that the two components of velocity are fixed on the lower boundary of this layer, so that the movement of lighter material dies out where the layer becomes very thin (rim synclines). Consequently, the new developing diapirs are mainly fed by material whose source lies in the nearly-undisturbed outer regions of the lower layer. If the boundary conditions on the upper and lower
~
h2
J
________________________________________________ ~~-‘
“~1” —
t~ 1.0
_~
-
1.14
-
1.46
Fig. 6. As for Fig. 5, except for the boundary conditions on the upper and lower boundary (see text for details),
.
the restrictions of the last model are no longer available. In this case, the distance between the different diapiric structuresis in good agreement with the results of the analytical solution, and the shape of the interface instability develops in a way that is dramatically different from that of its counterpart with the inverse boundary conditions. Another interesting phenomenon can be noted from Fig. 5. Although the mother diapir is already in a fully developed diapiric stage, the rise of lighter material can be re-activated. This can be seen in the local thickening of the thin pipe of fluid which connects the bottom layer with the upper part of the mother diapir, in agreement with observations of model experiments with different oils (Heye, 1978). Viscosities, layer thicknesses and time can be calculated for every model by dimensional analysis. The gravitational force is the same as every in nature, 3 in case.and Thethe visdensity/22contrast 0.1 the g cm cosity (poises)isand layer thickness h 2 (in cm) can be chosen independently. Then, time t (s) is given by t = t~ii2/h2.t can be directly taken from the plots. However, because of the very slow movement of the fluids, the Reynolds number (Re) has to be very small, and the restriction Re jh~/ii~I<< 1 between h2 and /22 can be derived.
4. Conclusions The two-dimensional Finite Element model presented here is able to follow the large deformations found in structures which originate from density instabilities. Calculations applied to single diapiric structures show the following results. (1) The shape of the growing interface instability is very dependent on the sign of the viscosity contrast and the ratio of thickness of the two layers. (2) For viscosity ratios clearly less than unity, an ~tinl growth can be found in the non(3) God agreement has always been found between analytical solutions and Finite Element calculations in the geometric linear stage of the interface mstabihty.
180
Results for calculations applied to diapiric structures, which are triggered by an initially isolated perturbation of the interface, show the following. (1) The distance between the different generations of diapirs is mainly controlled by the fastest-growing mode, as already predicted by the results of the linear analysis. (2) Migration of diapirs, renewed activation of upward movement, and large variations in the distance between different diapir generations, are phenomena associated in the geometric nonlinear stage of the interface instability, Acknowledgements The work was supported by the Deutsche Forschungsgemeinschaft.
References Danes, Z.F., 1964. Mathematical formulation of salt dynamics. Geophysics, 29: 414—424.
Heye, D., 1978. Experimente mit viskosen Flussigkeiten zur Nachahmung von Salzstrukturen. Geol. Jahrb., Reihe E, 12: 31—51. Hunsche, U., 1977. Modellrechnungen zur Entstehung von Salzdomfaniilien. Dissertation, Technische Universität Braunschweig, 102 pp. Nash, N.E. and Dabousi, O.B., 1978. Fluid dynamics model for salt-evolution. Tectonophysics, 47: 85—107. Prenter, P., 1975. Splines and Variational Methods. Wiley, New York, NY, 323 pp. Ramberg, H., 1967. Gravity, Deformation and the Earth’s Crust. Academic Press, New York, NY, 200 pp. Ramberg, H., 1968. Instability of layered systems in the field of gravity, land II. Phys. Earth Planet. Inter., 1: 427— 474. Schultz, M., 1973. Spline Analysis. Prentice-Hall, Englewood Cliffs, NJ, 156 pp. Selig, F., 1965. A theoretical prediction of salt dome patterns. Geophysics, 30: 633—643. Selig, F. and Wermund, E.G., 1966. Families of salt domes in the Gulf Coastal Province. Geophysics, 31: 726—740. Whitehead, J.A. and Luther, D.S., 1975. Dynamics of laboratory diapirs and plume models. J. Geophys. Res., 80:
705—717. Woidt, W.-D., 1979. Numerische Modellrechnungen zur Dynamic von Salzstockfamllien. Ph.D. Thesis (in preparation).
Physics of the Earth and Planetary Interiors, 21(1980)176—180 © Elsevier Scientific Publishing Company, Amsterdam — Printed in The Netherlands
FINITE ELEMENT MODELS OF DENSITY INSTABILITIES BY MEANS OF BICUBIC SPLINE INTERPOLATION W.-D. WOIDT
1
and HJ.
2
1 Inst itut für Geophysik, TU Braunschweig, Brunswick (Federal Republic of Germany) 2 Institu t für Geophysik, Goethe Universitàt, 6000 Frankfurt (Federal Republic of Germany)
(Received February 14, 1979; accepted for publication May 12, 1979)
Woidt, W.-D. and Neugebauer, H.J., 1980. Finite Element models of density instabilities by means of bicubic spline interpolation. Phys. Earth Planet. Inter., 21: 176—180. Density instabilities are often associated with large deformations in rocks, and result in diapiric structures of different shape and dimension. The Finite Element model presented here is able to follow the large deformations in a simple two-layer model with inverted density stratification. The developing interface instability gives interesting insights into the dependence of the rising diapiric structures on the model geometry, and the viscosity contrast between the two fluids. The development of single diapiric structures, as well as families of diapirs, are studied with the Finite Element model.
1. Introduction A model consisting of two superimposed viscous fluids with inverted density stratification, which develops spontaneously into a series of ridges and domes, has often been used to explain the dynamical behaviour of geological formations like salt domes, salt dome families, gneiss domes and, on a larger scale, continental rift systems and plumes. Such an unstable situation between fluids of different density and viscosity is generally called a Rayleigh—Taylor instability. As the viscosities of the fluids are very high if applied to geological timescales, the convective terms in the Navier—Stokes equation can be neglected, and nonlinearity only appears in the boundary conditions at the interface. For geophysical considerations, the temporal and spatial development of the interface between the two different fluids is of specific interest. Its shape can be studied as a function of the viscosity contrast and the boundary conditions at the interface. Figure 1 gives a short outline of the model geometry and the boundary conditions,
Analytical solutions of this two-dimensional model have been intensively studied in connection with the analysis of salt domes (Danes, 1964; Selig, 1965; Ramberg, 1968) and salt dome families (Selig and Wermund, 1966; Hunsche, 1977). These analytical solutions linearize the nonlinear term of the interface equation, and are therefore limited to infinitesimal disturbances. The linear model predicts an exponential growth of every sinusoidal disturbance in the interface. Disturbances of one specific wavenumber are found to show maximum growth rate, and therefore will determine the spacing of the developing interface structures. Because the linear model is limited to infmitesimal disturbances, it does not allow any assumption about the shape of the growing instability. The results of the analytical analysis can be applied to any artifical initial perturbation of the interface by decomposition into Fourier components. Every fundamental mode develops corresponding to its specific growth rate. After a time-increment, which is limited by the applicability of the linear analysis, the superposition of all independently
177
boundary conditions: z=h
1
:
0 or vI
=
:
Z=~(xt)
zs—h2 :
[ngn]~~
0
,
=
0
0. a’
=
=
0
0
~
=
•
0
0
Fig. 1. Two-dimensional model geometry and the boundary conditions.
growing components shows a wave-front propagation triggered by the initial perturbation. Such models have been successfully used to explain the phenomena of salt dome families, Starting from the results of the linear model, Nasir and Dabousi (1978) constructed a general solution by using’an infinite series of eigenfunctions. The infinite series solution is truncated to obtain the constants of integration from the boundary conditions along the interface, including the nonlinear term. Therefore their model makes it possible to study the shape of the growing instability. The results show the transition from the linear to the nonlinear development of the interface instability, but no diapiric stage. Up to now, diapiric structures have been studied mainly by model experiments using oils and glycerine (e.g., Whitehead and Luther, 1975) or with pasty materials in a centrifuge (e.g., Ramberg, 1967).
2. Numerical models In numerical calculations with the Finite Element method, the authors have determined that the main difficulty lies in following the large deformations in the two fluids. For this purpose, a deformable mesh, which follows the movement of the fluids, has only limited flexibility. The advantage of such a deformable mesh, however, is that the boundary conditions at the interface between different materials can be treated as natural, and therefore are easy to fulfil. The Finite Element model used here is based on the stream function approach, in which the condi-
tion of incompressibility is automatically implemented at the cost of fourth-order derivatives in the differential equation (Woidt, 1979). The interface between the two fluids is represented by discrete points, moving through a mesh, which remains stationary in space. Consequently, the abrupt change of the viscosities at the interface can only be described approximately in this model. Rectangular elements with high-order interpolation functions are used, which have to be of C’-continuity at least. Bicubic Hermitian polynomials (Prenter, 1975) as interpolation functions shows a very good order of accuracy, but unfortunately result in a large overall matrix system due to the sixteen unknowns per element. To follow the interface for as many time-steps as are necessary to reach its fully developed diapiric stage, within the capability of present day computers, requires basically small dimensions in the overall matrix. The number of unknowns can be considerably reduced by using bicubic spUme interpolation functions,which have a local support over sixteen elements and which are continuously differentiable twice (Schultz, 1973). The only unknown in this approach is the stream function itself. A typical mesh size of 16 X 25 rectangular elements results in 1768 unknowns and a semiband width of 72, using the bicubic Hermitian polynomials as interpolation functions in the Finite Element calculations. In comparison, the more complicated numerical algorithm, which makes use of bicubic spline interpolation functions, results in an overall matrix with 360 unknowns, and a semiband width of 49 for the same mesh.
178
The time integration scheme by which the interface is moved into its new position is a fourth-order Runge—Kutta or Adam—Bashford process. 40—100 time-steps are necessary to reach the fully developed diapiric stage. Lower-order schemes give wrong results because of the exponential or even over-exponential growth of the disturbance in the beginning of the dome’s development.
I
t~ 22.6
.~ “
20 ,
h1 18
~-.-
3. Results and discussion The first two plots (Figs. 2 and 3) show diapiric structures, which evolve from an initially sinusoidal perturbation of the interface. The only difference in the parameters of both models is in the ratio of viscosities, which is 100 : 1 in Fig. 2 and 1 :100 in Fig. 3. As it can be seen in the plots, the shape of the growing interface instability is very much dependent on the ratio of viscosities. Following the top of the rising diapir in Fig. 3, an over-exponential growth can be derived during the nonlinear stage, which is characteristic for viscosity ratios clearly less than unity. The thickness of the heavier upper layer is seven times that of the underlying lighter fluid. Because of this, the rising instability can develop its different stages without being influenced too early by the upper boundary. Such an early influence forces the upwelling lighter fluid to flow aside, forming a mushroom-shaped structure. The strong influence of this boundary on the development of the t~ ~ h
_________
.~—
8 _________
P
‘~
I Fig. 3. As for Fig. 2, except that the ratio of viscosities ~I/M2 = 0.01, and A = 2.96 h2. ,
instability can be studied in Fig. 4, where the thickness of both layers is about equal, and the ratio ofviscosities is 100 : 1. Compared with Fig. 2, the stage where the rising fluid forms a spherical bulb in Fig. 4 is suppressed by the upper boundary condition. In each case the wavenumber of the initial disturbance shows maximum growth as predicted by the linear analytical solution. Good agreement between the analytical solution and this numerical model has always been found for the linear stage of the interface instability. Figures 5 and 6 show diapiric structures where the initial triangular-shaped disturbance triggers the evolution of new diapirs From model experiments with oils it is known that viscosity ratios near to unity form the most fascinating instabilities (Heye,
~
1 h2:
17
13
~2
Fig 2 Model diapir which develops from an initially sinu soidal perturbation h1/h2 ~ = 100 and 3.A Free-slip = 13 h2 conditions on the boundary fixed (zero velocThe densityhold contrast forupper all models is 0.1and g cm ities) conditions on the lower surface. Time can be derived from the formula t = t*M2/h 2 (P2 in poises, h2 (cm) and t (s)). can be taken directly from the plots (see text).
1978). Therefore, a viscosity ratio of unity is chosen hi h2
__
________________________________________ .
.
—
—
.
—
Fig. 4. As for Fig. 2, except that the ratio of thicknesses h1/h2 = 0.92, and therefore A 5.36 h2.
179
boundary are exchanged, as has been done in Fig. 6, the dominant wavelength becomes much greater, and
~‘
h
2 1.0- 1.2- 1.32
Fig. 5. Model diapirs which are triggered by a triangular initial perturbation (width 6h2). Model parameters are ~ /~2= 1, hi/h2 = 7, and the length of the modelL = 48h2.
for these numerical calculations. The two examples differ only in the conditions on the upper and lower boundaries, which are free-slip on the upper boundary and non-slip on the lower boundary in Fig. 5, and vice versa in Fig. 6. The wavelength with maximum growth rate, as predicted by the linear solution, is expected, to dominate the distance between the different generations of diapirs,which is 3.7 times the thickness of the thin layer for Fig. 5. This theoretical prediction is in good agreement with the numerical calculations for the first generation of diapirs, far first as the early nonlinear stage is concerned. Later,asthe generation begins to grow into the direction of the mother diapir, and the distance between them becomes shorter. Finally, material streams into the second rim syncline, thereby distorting the first generation diapirs into a direction away from the mother diapir. It is interesting to note that the distance between first and second generation diapirs is much greater than that theoretically expected. This can be understood if one considers that the lighter fluid layer is too thin to feed a new diapiric structure every 3.7 h2, and that the two components of velocity are fixed on the lower boundary of this layer, so that the movement of lighter material dies out where the layer becomes very thin (rim synclines). Consequently, the new developing diapirs are mainly fed by material whose source lies in the nearly-undisturbed outer regions of the lower layer. If the boundary conditions on the upper and lower
~
h2
J
________________________________________________ ~~-‘
“~1” —
t~ 1.0
_~
-
1.14
-
1.46
Fig. 6. As for Fig. 5, except for the boundary conditions on the upper and lower boundary (see text for details),
.
the restrictions of the last model are no longer available. In this case, the distance between the different diapiric structuresis in good agreement with the results of the analytical solution, and the shape of the interface instability develops in a way that is dramatically different from that of its counterpart with the inverse boundary conditions. Another interesting phenomenon can be noted from Fig. 5. Although the mother diapir is already in a fully developed diapiric stage, the rise of lighter material can be re-activated. This can be seen in the local thickening of the thin pipe of fluid which connects the bottom layer with the upper part of the mother diapir, in agreement with observations of model experiments with different oils (Heye, 1978). Viscosities, layer thicknesses and time can be calculated for every model by dimensional analysis. The gravitational force is the same as every in nature, 3 in case.and Thethe visdensity/22contrast 0.1 the g cm cosity (poises)isand layer thickness h 2 (in cm) can be chosen independently. Then, time t (s) is given by t = t~ii2/h2.t can be directly taken from the plots. However, because of the very slow movement of the fluids, the Reynolds number (Re) has to be very small, and the restriction Re jh~/ii~I<< 1 between h2 and /22 can be derived.
4. Conclusions The two-dimensional Finite Element model presented here is able to follow the large deformations found in structures which originate from density instabilities. Calculations applied to single diapiric structures show the following results. (1) The shape of the growing interface instability is very dependent on the sign of the viscosity contrast and the ratio of thickness of the two layers. (2) For viscosity ratios clearly less than unity, an ~tinl growth can be found in the non(3) God agreement has always been found between analytical solutions and Finite Element calculations in the geometric linear stage of the interface mstabihty.
180
Results for calculations applied to diapiric structures, which are triggered by an initially isolated perturbation of the interface, show the following. (1) The distance between the different generations of diapirs is mainly controlled by the fastest-growing mode, as already predicted by the results of the linear analysis. (2) Migration of diapirs, renewed activation of upward movement, and large variations in the distance between different diapir generations, are phenomena associated in the geometric nonlinear stage of the interface instability, Acknowledgements The work was supported by the Deutsche Forschungsgemeinschaft.
References Danes, Z.F., 1964. Mathematical formulation of salt dynamics. Geophysics, 29: 414—424.
Heye, D., 1978. Experimente mit viskosen Flussigkeiten zur Nachahmung von Salzstrukturen. Geol. Jahrb., Reihe E, 12: 31—51. Hunsche, U., 1977. Modellrechnungen zur Entstehung von Salzdomfaniilien. Dissertation, Technische Universität Braunschweig, 102 pp. Nash, N.E. and Dabousi, O.B., 1978. Fluid dynamics model for salt-evolution. Tectonophysics, 47: 85—107. Prenter, P., 1975. Splines and Variational Methods. Wiley, New York, NY, 323 pp. Ramberg, H., 1967. Gravity, Deformation and the Earth’s Crust. Academic Press, New York, NY, 200 pp. Ramberg, H., 1968. Instability of layered systems in the field of gravity, land II. Phys. Earth Planet. Inter., 1: 427— 474. Schultz, M., 1973. Spline Analysis. Prentice-Hall, Englewood Cliffs, NJ, 156 pp. Selig, F., 1965. A theoretical prediction of salt dome patterns. Geophysics, 30: 633—643. Selig, F. and Wermund, E.G., 1966. Families of salt domes in the Gulf Coastal Province. Geophysics, 31: 726—740. Whitehead, J.A. and Luther, D.S., 1975. Dynamics of laboratory diapirs and plume models. J. Geophys. Res., 80:
705—717. Woidt, W.-D., 1979. Numerische Modellrechnungen zur Dynamic von Salzstockfamllien. Ph.D. Thesis (in preparation).