Summary
We present a numerical approach for solving 2-D mantle flow problems where the chemical composition changes abruptly across intermediate boundaries. The method combines a Galerkin-spline technique with a method of integration over regions bounded by advected interfaces to represent discontinuous variations of material parameters. It allows direct approximation of a natural free surface position, instead of a posteriori calculation of topography from the normal stress at the upper free-slip boundary. We formulate a model where a viscous incompressible fluid filling a square box is divided into layers (not necessarily horizontal) by advected boundaries, across which the density and viscosity change discontinuously. No-slip or free-slip conditions are assumed at the model sides. The suggested approach, being Eulerian, avoids the difficulties due to material discontinuities at intermediate boundaries, like the Moho or the Earth's surface, and is also free from the deficiencies of the Lagrangian approach, always resulting in mesh distortion. We present two geophysical cases analysed by this technique. The first case concerns the formation of sedimentary basins under the effects of heavy bodies sinking in the asthenosphere and of load due to sedimentary infills. The second case demonstrates the evolution of salt diapirs and shows how their growth is affected by a laterally inhomogeneous sedimentary layer. This numerical approach is well suited for problems of gravitational instability with discontinuities of density and viscosity across advected boundaries.
Introduction
The problems of gravitational instability involving distinct chemical layers are challenging in geophysical fluid dynamics. Motions of material interfaces separating geomaterials of differing material properties are essential to layered mantle convection, subduction of lithospheric slabs, ascent of mantle plumes, sinking of heavy bodies in the asthenosphere, salt diapirism and many other processes (Jacoby 1970; Römer & Neugebauer 1991; Ribe & Christensen 1994; Schubert, Anderson & Goldman 1995; Naimark & Ismail-Zadeh 1995). The advection of material interfaces was studied analytically for cases of small perturbations and displacements by Biot & Ode (1965), Chandrasekhar (1968), Ramberg (1968), Naimark & Yanovskaya (1976) and Naimark & Ismail-Zadeh (1994). A numerical approach is needed to examine finite displacements of material boundaries. Numerical schemes usually produce errors originating from discontinuities of physical properties, where a step function needs to be advected (Lenardic & Kaula 1993). To minimize such errors, methods of representing these discontinuities were developed by Christensen (1982, 1992) and Naimark & Ismail-Zadeh (1995, 1996).
The problem is simple for the case of discontinuous density, because it enters the right-hand side of the momentum equation. The method of tracer chains suits the purpose, because it reduces the problem of density discontinuities to the integration of relevant terms along the curves (tracer chains).
The problem is much more difficult for the case of discontinuous viscosity. It enters the momentum equation and hence, in the Eulerian approach, must be locally smoothed over several grid steps (e.g. Woidt 1978; Christensen & Yuen 1984; Christensen 1992; Naimark & Ismail-Zadeh 1995, 1996). Additionally, a global filter technique (Lenardic & Kaula 1993) is used to remove overshoots and undershoots of an advected step function. Dense enough grids are needed to approximate discontinuities by smooth functions, which requires much computer resources.
We present a numerical approach for solving 2-D Stokes' flow problems where physical properties (density and viscosity) change discontinuously across advected boundaries. The approach combines a Galerkin-spline technique with a method of integration over advected layers, where a finite-dimensional space of spline weights is used together with Cartesian coordinate representations of discontinuous viscosity terms. It allows approximation of the natural shape of a free surface, instead of a posteriori calculation of its topography from the normal stress at the upper free-slip boundary. We present two geophysical cases analysed by this technique. The first case concerns the formation of sedimentary basins under the effects of heavy bodies sinking in the asthenosphere and of loads due to sedimentary infills. The second case demonstrates the evolution of salt diapirs with laterally homogeneous and inhomogeneous overburdens of sediments.
Model Concept
Fig. 1 illustrates the rectangular model region Ω: 
and H are model width and depth; a Newtonian fluid with variable density ρ and viscosity μ fills this region. Curves ℒe= 1,2,…, E divide the model region Ω into several subregions Ωe, e= 1,2,…, E + 1. We assume that each curve ℒe is closed or starts and terminates at the boundary of Ω and has no self-intersections; however, different curves can intersect each other. Fig. 1 shows two curves, ℒ1 and ℒ2, and three subregions Ω1, Ω2 and Ω3. In what follows, we consider one curve ℒ for simplicity, though the number of curves can be arbitrary. We also use dimensionless forms of equations governing the model, so that after the appropriate change of variables the model region Ω occupies the square 0 ≤x≤ 1, 0 ≤z≤ 1.
Geometry of the model for the case of two interfaces.
Introduce the following notation: 
where ψ(x, z, t), ψ(x, z, t) and A(x, z, t) are functions having continuous derivatives entering in the notation.
We seek the stream function ψ, density p(x, z, t), viscosity μ(x, z, t) and the family of curves ℒ: x=x(q, t), z = z(q, t) (q is a parameter of points on a curve, 0 ≤q ≤ Q) satisfying the differential equations (g is the acceleration due to gravity)
the impenetrability and free-slip boundary conditions
and initial conditions at t=t0
The first equation is the 2-D Stokes equation represented in terms of the stream function ψ, the second and third equations describe the transfer of density and viscosity with the flow, and the remaining equations determine trajectories of points x(q, t) and z(q, t) located at t0= 0 on the curve ℒ0.
We define a weak solution of the problem, that is, a solution satisfying an integral relation rather than the equation itself. Let us multiply the first equation in (1) by a function ψ(x, z, t) satisfying the same boundary conditions as ψ(x, z, t), integrate by parts the left- and right-hand sides of the product twice and once, respectively, and observe that the integral over the model boundary vanishes. Multiply the second and third equations in (1) by functions ϑ and ζ, respectively, and integrate the results. A weak solution of the problem stated above is the set of functions ψ(x, z, t), p(x, z, t), μ(x, z, t), x(q, t) and z(q, t) satisfying the above boundary and initial conditions and the following equations:
where ψ(x, z) ε, ϑ(x, z) εℜ and ζ(x, z) εℜ are any functions (called test functions) from sets and ℜ properly chosen.
Numerical solutions are obtained in the form of weighted sums of basic bicubic splines. However, bicubic splines, being excellent for the case of smooth unknown functions, become inadequate when these functions are discontinuous. To preserve the accuracy of spline representations for cases of discontinuous unknowns, we suggest the method described below.
Let us represent unknown functions ρ(x, z, t) and μ(x, z, t) as sums of two functions, one smooth and the other constant over Ω1 and Ω2:
where ρ1(x, z, t) and μ1(x, z, t) have continuous first and second derivatives, whereas ρ0(x, z, t) and μ0(x, z, t) take on constant values in Ω1, and Ω2:
where 
, 
, 
and 
are functions of time, but do not depend on x and z. Let us substitute representations (3) into the first relation in (2) and obtain the result
and in the interior of any region Ω, (see Fig. 1)
because 
in this interior. These equations, together with
and with boundary and initial conditions described above, define a weak solution for the case of discontinuous density and viscosity.
Numerical Method
Approximations to unknown functions ψ, ρ1 and μ1 are represented as linear combinations of basic bicubic splines with unknown coefficients (here and below we assume summation over repeated subscripts taking on the following values: i, k, m= 1, …, I; and j, l, n= 1, …, J):
where spij(x, z) and 
are basic bicubic splines satisfying the required boundary conditions. These splines are constructed from basic cubic splines: spij(x, z) = si(x)sj(z) and 
The basic cubic and bicubic splines used here were described by Naimark & Malevsky (1986). The curve ℒ is approximated by a polygon whose vertices have coordinates xβ(t), zβ(t), β= 1, …, B. These vertices are located on ℒ0 at t = t0. Let us substitute the above representations into eqs (4) and (5) and integrate forms involving products of basic splines and their derivatives. This results in sets of linear algebraic equations for unknowns ψij, ordinary differential equations for 
and 
:
Coefficients Cijkl are sums of three terms, 
, where the first term is obtained from μ1 by substituting its spline representation into the first integral in (4), rearranging sums, and integrating products of splines and their derivatives. The result takes the form
where
Here (…)(p) denotes the derivative of order p of a function (…) and the zero-order derivative is the function itself. The terms 
and 
are obtained by integrating products of splines and their derivatives over regions Ω1 and Ω2, which results in the forms
We see that elements 
depend on the continuous term μ1, but are independent of the curve ℒ. On the other hand, elements 
and 
depend on the curve ℒ and on the constants 
and 
but are independent of the continuous term μ1.
Coefficients Fijkl on the right-hand side of the first equation in (7) are obtained by integration:
The term Ψkl is obtained from the last two integrals on the right-hand side of eq. (4), where θ is set to sk(x)si,(z). The sum of these integrals takes the form
explained in detail by Naimark & Ismail-Zadeh (1995).
Coefficients Gijkl and Eijkl entering the second and third equations in (7) are also calculated by integrating basic splines and their derivatives:
where 
and 
are obtained from 
and 
in eqs (9) with si(x), sk(x), 
, Sj(z), si(z) and 
replaced by 
, 
, sm(x), 
, 
and sn(z), respectively.
The unknowns to be found from eqs (7) are the following: ρij(ts), μij(ts), ψij(ts), xβ(ts) and zβ(ts), s = 1, 2, …, S. The second. third, fourth and fifth relationships in (7) constitute the set of ordinary differential equations (ODE) for unknowns ρij, μij, xβ and zβ. We solve this set by the fourth-order Runge—Kutta method. The right-hand sides of these equations include unknowns ψij found from the first set of equations in (7). Initial values pij(t0) and μij(t0) are derived from the conditions
by using spline interpolation programs. Let us describe the calculation of the right-hand sides. Assume that the unknowns have been calculated at t = t5. Use eqs (8)–(10) to find the matrix Cijkl and eqs (11) and (12) to compute the right-hand sides of the first set in (7). Solve this set for φij and use the values so found, together with eqs (13), to calculate the right-hand sides of the above ODE.
Coefficients (9), (11) and (13) can be computed once and used in all calculations. Certain difficulties arise in eqs (10). The integrals in eqs (10) depend on the curve ℒ, changing with time. We reduce calculations of the forms (10) to direct integration of polynomials over regions bounded by the curve ℒ and model boundaries; these polynomials are products of splines and their derivatives. To clarify this reduction, consider Fig. 2. This figure shows part of a rectangular grid, the curve ℒ passing through it, and parts of Ω1 and Ω2. The integrals in eqs (10) can be treated as sums of integrals over all grid squares. They are easily computed for squares that are not intersected by ℒ. When a grid square is intersected by the curve, as shown in the figure, the integration is carried out over a region bounded by the square sides and the part of ℒ in this square. This is done by summing up integrals over all trapezoids bounded by the edges of the polygonal curve ℒ, vertical lines and horizontal segments of the lower grid square side; one such trapezoid is shaded in Fig. 2. The integral over each trapezoid is computed directly by repeated integration of polynomials. Obviously, horizontal trapezoids are treated instead of vertical ones when ℒ passes through a square from top to bottom. The crucial part of this procedure is the analytical integration over regions whose boundary includes ℒ, which can intersect grid squares in many ways.
A sketch of the curve passing through a rectangular grid. A current trapezoid T is shaded. The integrals in eqs (10) are calculated as sums of integrals over all trapezoids in the grid square intersected by the curve.
Let us discuss the conditions at advected material boundaries. Physically, velocity and stress are continuous across these boundaries. It follows that the viscosity discontinuities across the same boundaries lead to discontinuity of the second invariant of strain rate. In the approach described here the second derivatives of the stream function ψ are continuous, which can look inconsistent with the physical conditions. However, the suggested algorithm leads to fitting these conditions with a continuous ψ, which results in a ψ smoothed over a chosen grid. Overshoots and undershoots of a stream function so found are not as great as might be expected, because the algorithm fits these conditions for second partial derivatives of ψ, rather than for ψ itself.
It is, of course, impossible to formulate strictly conditions at free boundaries in the stream function approach, because the order of equations changes when viscosity equals zero. However, when viscosity at one side of a surface is sufficiently low compared to that at the other side, the algorithm still works and leads to correct results. This is verified here in the case of isostatic adjustment in a layered medium. We call boundaries of this kind ‘free’, to emphasize that they are approximations to physical free boundaries.
Let us note that an interface ℒ, being advected, stretches or compresses, so that the distance between adjacent vertices of its polygonal representation can become too long, and computations lead to erroneous results or deteriorate. To avoid this, we periodically update the polygonal line ℒ. Denote by |ℒ| the length of ℒ and by n the number of its vertices. New vertices are placed on ℒ with the step 
along it. This results in an almost uniform spacing of vertices on ℒ at any time step.
Verification Of The Method
An exact solution of eqs (4)–(6) is unknown even for the simplest cases and boundary conditions. Previous methods were tested for eqs (4) and (5) separately (Ismail-Zadeh, Naimark & Lobkovsky 1996; Naimark & Ismail-Zadeh 1996). The suggested algorithm and codes were verified by comparing numerical results with experimental data and analytical results from the linear theory of the Rayleigh—Taylor instability.
We studied a model of viscous layers with stable density stratification and used data employed in the experimental isostatic test performed by Ramberg (1968). In Ramberg's physical experiment a heavy syrup was supporting a less dense layer of a silicone putty. Initially, a sinusoidal deflection was prescribed at the free surface, whereas the interface between the substratum of syrup and silicone putty was straight and horizontal. The experiment showed that the free surface tended to flatten, and the interface was deflected into a wave. After passing a maximum amplitude, the secondary wave at the interface began to flatten until a stable equilibrium was attained.
We modelled this experiment as follows. The model region was divided into three layers by two curves: 
, and 
. The densities, viscosities and thicknesses of the layers are shown in Table 1. Curve ℒ1 approximates a free surface; the wavelength of the perturbation λ equals model width L. We used two rectangular grids (20 × 25 and 46 × 48) and obtained the same results. The pattern of layers obtained numerically is very close to the experimental results of Ramburg (1968).
Nomenclature of and values used in the model.
Fig. 3 shows an amplitude of wave versus time at the upper boundary and the interface obtained from measurements and the numerical test. The solid lines calculated by the suggested method are indistinguishable from theoretical curves predicted by the linear theory of gravitational instability (Ramberg 1968; Naimark & Yanovskaya 1976). We also see that the computed curves agree closely with experimental data.
Amplitude versus time for a model of isostatic adjustment in a layered medium. Solid lines show results obtained by the suggested numerical model. Points and crosses represent results from three experimental runs (Ramberg 1968, Fig. 25).
Efficiency Of The Method
We analysed two numerical approaches, one proposed previously by Naimark & Ismail-Zadeh (1995), hereinafter called ‘old’, and that suggested in this paper, called ‘new’, for the case of a model consisting of three layers. This model is sketched in Fig. 4(a) with the grid used in the calculations. Region 1 has zero density and low viscosity; region 2 is highly viscous, heavy and thin. We call this region thin because it contains two or three grid levels in the z direction. Region 3 is less dense and viscous than region 2. The values of the dimensionless model parameters are presented in Fig. 4(a). Figs 4(b) and (c) show the positions of layer boundaries calculated by the ‘old’ and ‘new’ methods at different times.
Testing the ‘old’ and ‘new’ numerical approaches for the case of a layered model with a heavy, highly viscous layer and a ‘free’ surface. (a) A sketch of the model and the grid chosen for calculations. The curves depict initial positions of boundaries between layers. (b) The positions of upper and lower boundaries, shown separately and appropriately scaled, at times 479.5 (‘old’, solid lines) and 473.7 (‘new’, dashed lines). (c) The same for times 2525.0 (‘old’) and 2421.0 (‘new’).
We see that the ‘old’ method leads to erroneous perturbations of advected boundaries. It is seen that these perturbations grow with time, while they are absent when the positions of the boundaries are calculated by the ‘new’ method. Subsequent calculations by the ‘old’ method lead to an almost instant deterioration of the pattern when time reaches a threshold (in this case about 3000), whereas the layered pattern remains adequate with the ‘new’ method.
The deficiency of the ‘old’ method in cases of thin layers can be explained as follows. When the number of grid points across a layer is small, overshoots and undershoots of viscosity lead to large errors in computing the stream function, hence to erroneous velocities controlling advection of boundaries. The values of viscosity at some points can even become negative (as in this test where the viscosity changes from 1 to 1000 across the upper boundary). Naturally, negative viscosity results in erroneous velocities. This effect leads to the deterioration of the overall pattern. On the other hand, there are no overshoots and undershoots of viscosity in the ‘new’ method. Viscosity remains constant in each of the advected regions.
It is possible to obtain correct results by the ‘old’ method with denser grids, advanced approaches to smoothing viscosity across the boundaries, and very small time steps. However, the computer resources required will become much greater than those needed in the ‘new’ method.
Sample Calculations
Sinking of heavy bodies and sedimentary basin formation
In the magmatism—eclogitization mechanism of sedimentary basin evolution (Lobkovsky . 1993; Ismail-Zadeh . 1996) eclogitic bodies evolve from magmatic melts accumulated at a depth of 60–80 km in a post-rift phase. These bodies, being denser than the surrounding material, sink in the astheno-sphere and induce viscous flows that change surface topography and lead to the formation of sedimentary basins. At least three boundaries where physical properties are discontinuous should be introduced in this model: ℒ1, the ‘free’ surface; ℒ2, the sediment/basement interface; and ℒ3, the heavy body/asthenosphere interface. Naimark & Ismail-Zadeh (1995) described a similar model where sediments were absent and where the surface topography was calculated a posteriori from the normal stress at the upper free-slip surface. The present model includes sedimentary infill and viscosity discontinuities.
The model of sedimentation, chosen here for its simplicity, is based on two assumptions: (1) the depression is filled by sediments instantaneously; (2) this filling stops when the upper level of sediments achieves the value prescribed by a constant a. Sediments in this model appear ‘from nowhere’; we do not consider processes leading to sedimentation. However, the law of conservation of mass is not violated: the initial mass remains constant, and the new mass is added by sedimentation processes. Note that it is quite easy to incorporate other models by introducing the rate of sedimentation or considering mechanisms based on the shape of the free surface.
According to this model, boundaries ℒ1 and ℒ2 are initially (at t= 0) coincident with straight horizontal lines z = b, where b < 1 is a constant. Boundary ℒ1 is deflected by the viscous flow, and boundary ℒ2 is calculated in the following manner. Denote by 
the function representing curve ℒ1 at time t. Find 
and 
at a current time t and put 
where a, 0 ≤, ≤ 1, is a constant controlling the filling of the basin with sediments. The curve ℒ2 is computed at any time t from the condition
Hereinafter all variables are dimensionless, unless otherwise stated. The following initial geometry was assumed: the ‘free’ surface z= 0.77 and the heavy ellipse centred at x= 0, z= 0.5 with vertical and horizontal semi-axes 0.03 and 0.4, respectively. The density above the ‘free’ surface was 0.0, within sediments 2.5, in the asthenosphere 3.5, and 4.0 in the ellipse. The viscosity was 1.0 above the ‘free’ surface, 100.0 in the asthenosphere and sediments, and 110.0 within the ellipse. Numerical tests showed that viscosity variation above the ‘free’ surface from 1.0 to 10−2 resulted in very small changes in numerical solutions. The constant a was taken equal to 0.7: at any time the depression is filled with sediments to 0.7 of its depth. Calculations were made with a rectangular 20 × 25 grid.
Fig. 5 shows four snapshots of a flow at different times. Each snapshot consists of the lower and upper panels. The lower panel depicts the position of the ‘free’ surface and heavy body. The upper panel shows the vicinity of the ‘free’ surface stretched in the vertical direction to make sedimentary infill and changes of topography clearly visible. We see how the ‘free’ surface, initially flat, deflects under the actions of the sinking heavy body and of loads due to sediments.
Four phases of sedimentary basin evolution under two effects: the flow produced by a sinking heavy body (shaded) and load due to a sedimentary infill (shaded). Two panels illustrate each phase. The lower panel depicts the position of the ‘free’ surface and heavy body. The upper panel shows the vicinity of the ‘free’ surface stretched in the vertical direction to make sedimentary infill and changes of topography clearly visible. We see how the ‘free’ surface, initially flat, deflects under the actions of the sinking heavy body and of loads due to sediments.
A test with a= 0 (no sediments) leads to similar patterns, but with less subsidence. Fig. 6 shows two subsidence curves: for the cases a= 0 (curve 1) and a= 0.7 (curve 2). The subsidence was calculated as a depth from the initial position of the ‘free’ surface to the deepest point of the deflected boundary. We see that sedimentary loads can increase the basement subsidence by a factor of 2 or more.
Modelled subsidence curves: in the absence of sedimentary loads (1) and with sediments (2).
Evolution of salt diapirs
Salt diapirism is another process involving viscous flows with material boundaries. Salt tectonics is quite important from the practical point of view, because various types of hydrocarbon traps are closely associated with salt domes (Talbot 1992). Numerical models of salt diapirism were extensively studied by Woidt (1978), Schmeling (1987), Römer & Neugebauer (1991), Poliakov & Podladchikov (1992), Zaleski & Julien (1992), Poliakov . (1993), Podladchikov, Talbot & Poliakov (1993), Keken . (1993) and Daudré & Cloetingh (1994).
Natural salt structures have various shapes (Jackson & Talbot 1986, Volozh, Groshev & Sinelnikov 1994), which strongly depend on the thickness of the salt layer and that of the surrounding overburden (Schmeling 1987) and on the horizontal gradient of loads due to sediments (Poliakov . 1993). We present two cases of a salt layer evolution: model A, with the ‘balloon on a string’ geometry, and model B, where nappes superimposed on the sedimentary overburden lead to asymmetric diapirism. The square model box is 15 km long and 5 km deep. This box is divided into 76 × 26 rectangular elements in the x and z directions, respectively.
Model A
A salt layer 0.5 km thick at the bottom of the model is covered by a sedimentary overburden 2.5 km thick. The salt/sediment interface is initially perturbed by a peak of cosine shape with amplitude 0.2 km and length 0.57 km. The viscosities and densities are 1020 Pa s and 2.3 × 103 kg m−3 for the overburden and 1018 Pa s and 2.2 × 103 kg m−3 for the salt.
Fig. 7 shows the evolution of a diapir evolved from the initial perturbation in 47 Myr. The shapes of salt structures closely agree with classical cases of the Rayleigh—Taylor instability with high viscosity contrasts and a thin lower layer.
Evolution of salt diapirs, model A. Viscosities and densities are as follows: 1020 Pa s and 2.3 × 103 kg m−1 (overburden, medium shading) and 1018 Pa s and 2.2 × 103 kgm−3 (salt, heavy shading). (a) t= 0; (b) t= 27.4 Ma; (c) t= 33.7 Ma; (d) t= 36.6 Ma; (e) t= 41.8 Ma; (f) t= 47.1 Ma. The timescale used is t*=μ*/(ρ*gH)= 33.3 yr, where μ*= 1017 Pa s and ρ*= 2.3 × 103 kg m−3. Flow velocities are shown by arrows. The velocity scale is given at the top of each figure.
Model B
This model presents salt motions in the presence of laterally asymmetrical loading. We feel that nappes of sediments can lead to asymmetrical shapes of salt structures (C. Talbot, personal communication, 1996). The nappe of sediments in model B has maximum thickness 0.99 km, viscosity 1020 Pa s and density 1.9 × 103 kg m−3. Fig. 8 shows the evolution of the resultant salt structure. The nappe of sediments was imposed on the overburden with the growing symmetrical diapir (Fig. 8a). The velocity of nappe sinking in the overburden is greater than the rate of diapiric growth. This is clearly seen from velocities presented in Figs 8(a) and (b). Fig. 8(b) demonstrates also how the shape of the diapir becomes slightly asymmetrical. When the nappe attains its equilibrium, the rate of diapiric penetration increases (Fig. 8c). Subsequent phases of diapiric evolution are shown in Figs 8(d) and (e). It is seen that the diapir remains only slightly asymmetric. However, even this minor asymmetry leads to a quite asymmetric shape of the diapir in its subsequent evolution (Fig. 8f).
Evolution of salt diapirs under the effect of laterally inhomogeneous sedimentary loads, model B. Three layers are present: salt (heavy shading), overburden (medium shading) and nappe of sediments (light shading). Viscosity and density of sediments' nappe are 1020 Pa s and 1.9 × 103 kg m−3. Viscosities and densities of other layers, the timescale and velocity representations are the same as in Fig. 7. (a) t= 0; (b) t= 0.3 Ma; (c) r= 3.3 Ma; (d) t= 6.4 Ma; (e) t= 10.2 Ma; (f) t= 30.5 Ma.
Discussion And Conclusions
The method suggested here results in advection of step functions (density and viscosity) free of overshoots and undershoots: the values of ρ0 and μ0 in regions bounded by interfaces remain unchanged by definition. However, the method has its limitations, because the difficulty of representing discontinuous changes of physical properties shows up elsewhere. When viscosity is discontinuous across an interface, natural boundary conditions (continuity of stress and velocity) result in discontinuous strain rate, expressed in terms of the second derivatives of the stream function. However, the stream function, being a spline in the assumed approach, must have continuous second derivatives. As a result, the Galerkin method yields the strain rate locally smoothed in the vicinity of the interface. The velocity in the vicinity of the interface is continuous but has a sharp variation. This leads to overshoots and undershoots of the velocity. However, these overshoots and undershoots are not great, because second derivatives of ψ rather than ψ itself are smoothed at the interface.
Numerical tests of the previous method (Naimark & Ismail-Zadeh 1995) show that errors in advected step functions grow with time. A smoothing technique can reduce these errors, but they always tend to increase; they can be treated as perturbations giving rise to new instabilities. Sometimes this leads to erroneous patterns that look like mixing and can result in wrong conclusions. In the suggested method, overshoots and undershoots of velocities do not grow with time; moreover, numerical tests show that these errors tend to decrease on attaining a certain level. From this viewpoint, the present method is more stable than other Eulerian numerical methods involving advection of step functions. Testing of the method for the case of a thin layer shows its advantages.
A numerical model of isostatic adjustment of a layered medium demonstrates a very close agreement with experimental data and results predicted by the linear theory of gravitational instability. The method allows one to take into account the appearance of additional structures bounded by material interfaces, such as nappes of sediments or sedimentary infills.
Acknowledgments
The authors are grateful to C. Talbot and Yu. Volozh for useful discussions on salt tectonics and to D. Yuen for discussing applications of the method. We would like to thank P. van Keken, A. Malevsky and H. Schmeling for valuable reviews of the initial version of the manuscript and for discussions. ATI is very grateful to R. Nicolich and L. Cernobori, DINMA, University of Trieste, and to E. Aureli and P. Öster, PDC/NADA, Royal Institute of Technology, Stockholm, for computing facilities used to perform a large part of the calculations. ATI was supported by the Swedish Institute and Stockholm University during his stay in PDC/NADA. This research was supported by INTAS (grant 94-1099) and RFBR (grants 96-05-64356, 97-05-65415).
References
