Numerical simulations of the propagation path and the arrest of £uid-¢lled fractures in the Earth

SUMMARY We use a boundary element method to study the growth and quasi-static propagation of £uid-¢lled fractures in regions with inhomogeneous and deviatoric stresses. The wholesale migration of fractures due to their opening at one end and closing at the other can be simulated when using a ¢nite £uid mass contained in a fracture and considering £uid compression or expansion with changing fracture volume; these fractures are driven by stress gradients and by the density di¡erences between the £uid and the surrounding rock. Contrary to commonly held beliefs, the fracture growth and the propagation directions are not controlled only by the direction of the principal stresses, but also by tectonic stress gradients, apparent buoyancy forces and the length of the fractures themselves. The models help to explain the formation of sills, the lateral migration of magmas under volcanoes and the absence of volcanoes under the shallow parts of the Nazca plate. of the maximal strain-energy release and is determined from the strain-energy release in the three trial directions and inverse parabolic interpolation.


INTRODUCTION
Fluids play an important role in the geochemical and geophysical structure and mechanics of the Earth's lithosphere. For instance, magmatic dykes are probably the most common major structural element of the oceanic crust. Dyke and sill emplacements at volcanoes produce signi¢cant stress ¢eld inhomogeneities and are the major volcano-forming structure. A more recent recognition is the fact that meteoric £uids can penetrate at least 10 km into the continental crust, and the horizontal scale of meteoric £uid transport is now known to extend hundreds of kilometers (e.g. Torgersen 1991). The meteoric £uids in the crust seem to have a strong in£uence on the stability of rocks and may play a role in the triggering of earthquakes (e.g. Hickman et al. 1995). Furthermore, work by petrologists over the last few decades has lead to an understanding and recognition of the importance of £uid £ow in metamorphism (e.g. Ferry 1994;Rumble 1994), i.e. at deep crustal levels. It is thus apparent that any description of lithospheric dynamics is incomplete without an understanding of the in situ £uid properties and their transport mechanisms.
The`traditional' models of £uid transport in the lithosphere are either diapiric ascent or £uid £ow through networks of interconnected static microfractures (porosity networks). The diapiric ascent model is mostly accepted for magma batches that rise slowly in the upper mantle and to explain plutonic structures found in the continental crust. However, the diapiric ascent model is also controversial since numerical modelling indicates that the ascent is so slow that solidi¢cation of the magma may occur below the depth where some plutons have formed (see Clemens & Mawer 1992), and since mantle xenoliths provide evidence that a fast transport mechanism must exist in some cases to bring mantle magma to the surface.
Flow through porosity networks is an accepted model, for example, for £uid transport in sediments at low con¢ning pressures. However, recent studies on geopressurized sediments in the Gulf of Mexico indicate that £ow through porosity networks cannot explain all aspects of the observed large-scale £uid £ow in the geopressured zone (e.g. Nunn 1996). Porosity network models may oversimplify the situation at high con-¢ning pressure if the interaction between £uid-¢lled fractures and the elastic rock is not considered.
The propagation of isolated £uid-¢lled fractures and magma dykes is an important third model for the mechanism of £uid transport through the lithosphere and possibly the asthenosphere. Fluid-¢lled fractures can propagate much faster than diapirs and may thus overcome some of the shortcomings of the diapiric ascent model. The idea of propagating fractures also gets support from geological observations, which provide evidence that fracturing is a concurrent process during diapiric magma ascent; for example, magma dyke intrusions measured at the Oman Ophiolite show that magma-induced fracturing is a continuous process and one that spans the evolving structure of a dynamic diapir at the Moho level (Nicolas et al. 1994). Additionally, the propagation of isolated £uid-¢lled fractures is of interest when trying to consider the coupling between fractures and rock in porosity network models. In this paper we numerically study the propagation of £uid-¢lled fractures in elastic rock. Weertman (1971a) was one of the ¢rst to analyse quantitatively the break-open, heal-shut, rise-up behaviour of a single £uid-¢lled fracture in the interior of a plate. This mechanism was postulated to be relevant for the sinking of water-¢lled glacier crevasses and for the rising of natural hydraulic fractures and magma-¢lled dykes in the lithosphere (e.g. Weertman 1971bWeertman , 1973Weertman , 1980Weertman & Chang 1977;Secor & Pollard 1975;Pollard 1976;Pollard & Mueller 1976), and has been demonstrated in laboratory experiments of £uid-¢lled fractures in gelatin (Takada 1990;Heimpel & Olson 1994). [The downward £ow of water through temperate glaciers seems to be dominated, however, by the melt-driven down-cutting of horizontally £owing water through crevasses (e.g. Fountain & Walder 1998).] Although the analytical models of buoyancydriven fracture propagation have elucidated how upward fracture propagation may occur, they are unable to quantify propagation paths of £uid-¢lled fractures in inhomogeneous media and stress ¢elds. Our numerical approach is well suited to simulate such situations and is therefore complementary to the early analytical models.
We extended the displacement discontinuity method of Crouch & Star¢eld (1983) to study £uid-¢lled fracture propagation. In an earlier paper (Dahm 1996) we described how we extended this boundary element method for crack problems in arbitrary layered media and with surface topography. Here, for simplicity, we will not describe the layer and surface topography problem; we will only concentrate on the major steps necessary to handle quasi-static £uid-¢lled fracture propagation in the lithosphere.
The numerical simulation of quasi-static fracture growth has been performed previously (e.g. Sempere & MacDonald 1986;De Bremaecker & Swenson 1990;De Bremaecker & Wei 1994;Wei & De Bremaecker 1995a,b;Shen & Stephansson 1995). However, in our terminology these authors simulated empty fracture growth and not wholesale fracture motion, that is, the displacement of a whole fracture from one position to another (denoted as fracture propagation in the following). They disregarded the possible volume expansion or compression that may occur when £uid-¢lled fractures propagate; however, this must be considered for the simulation of £uid-¢lled fracture propagation. We have added two simple but crucial concepts in our approach to make our method well suited for the study of £uid-¢lled fracture propagation. One is the introduction of inequality constraints for opening dislocations at the fracture, in order to avoid a possible interpenetration of fracture surfaces under compressive stress. The other is the linear incorporation of a £uid pressure change with changing fracture volume. As a consequence, we have added another parameter to the problem, the initial volume per unit length perpendicular to the fracture, which has an in£uence on the propagation paths and lengths of £uid-¢lled fractures. The fracture propagation direction is estimated using a maximal strain-energy release rate criterion.
In the ¢rst part of this paper we describe the numerical method, in the second part we discuss the principles of £uid-¢lled fracture propagation and in the third part we predict the e¡ect of surface topography and asthenospheric mantle £ow above a subducting slab on the propagation of £uid-¢lled fractures and dykes. In a separate paper we study the velocity of £uid-¢lled fracture propagation (Dahm 2000).

NUMER ICAL METHOD
Our method is based on the displacement discontinuity method of Crouch & Star¢eld (1983). As is typical for boundary element methods, the computation is divided into two steps: the ¢rst is to solve the boundary value problem for prescribed stresses or displacements on the fracture surface, ending with a set of dislocations occurring at the fractures; the second is to use these dislocations to calculate displacements and stresses at arbitrary points in the studied area. The ¢rst step is the most di¤cult one and is described brie£y here.
According to Crouch & Star¢eld (1983) the boundary value problem at a fracture can be written in compact linear form as where *u N j and *u S j are the unknown normal and shear dislocations at the fracture segment j, respectively, B N i and B S i are the normal and shear tractions prescribed at the midpoint of segment i, C kl ij are the corresponding in£uence coe¤cients or Green's functions between observation and source segment i and j, and M is the number of straight-line segments used to approximate the fracture. The coe¤cients C kl ij for a full and a half-space are given in Chapter 5 of Crouch & Star¢eld (1983). The prescribed traction B S i is determined from the shear stress p S at the fracture, while B N i is a superposition of the components determined from the normal stress p N and from the internal £uid pressure P. The exactly determined linear system of equations (1) can be solved using a standard method. We solve (1) under the additional constraints using the non-negative least-squares method of Lawson & Hanson (1974). This ensures that the interpenetration of opposite fracture surfaces is impossible even under compressional external stresses. An important modi¢cation of the method is to include the £uid compression or expansion with changing fracture length and volume, and to consider its e¡ect on the internal £uid pressure. The relative volume change is given by (A{A 0 )/A 0 , where A 0 and A are the cross-sectional areas of the fractures with lengths l and lzdl, respectively. For example, the relative volume increase of gas-free magma may be about 10 per cent for a dyke ascending 100 km from depth. The corresponding density decrease is about 9 per cent. Note that £uid compression or expansion with changing volume is unimportant for non-propagating fractures. Since we study the quasi-static growth of fractures, the pressure change *P in response to a relative volume change is acting with equal magnitude at every segment of the fracture. The pressure change at segment i can be interpreted as a superposition of components resulting from volume change at every segment j, where K is the bulk modulus of the £uid and a j is the halflength of segment j. We use this simple isothermal equation of state, and disregard possible phase transitions. Considering the e¡ect of £uid compression by adding (3) to the left side of prescribed tractions B N i in (1) gives The equation for B S i is unchanged. Eq. (4) is still linear. A di¤culty may arise from possible di¡erent magnitudes of K and B N i ; K is usually in the range of a few gigaPascal, while B N i may be much smaller. We therefore use double precision ¢elds in our program; this seems to be su¤cient for most practical problems. If, as for a gas-¢lled fracture, K is very small and approximately zero in the limiting case, we call it an empty fracture.
To simulate quasi-static fracture growth and fracture closure we incrementally increase or decrease the fracture length. The quasi-static formulation is appropriate for £uid-¢lled fracture propagation since propagation velocities are much smaller than the elastic wave velocities (in the range 0.2^2 m s {1 , Brandsdo¨ttir & Einarsson 1979 and unpublished data) so that dynamic e¡ects can be neglected. We use a maximum strainenergy release rate criterion (generalized Gri¤th criterion) to estimate whether and in which direction a fracture will grow. A similar criterion has been used by De Bremaecker & Wei (1994) and Wei & De Bremaecker (1995a,b) to simulate mixed-mode fracture growth under shear load and compression. Other authors have used a maximum tensile or maximum circumferential stress (e.g. Sempere & MacDonald 1986;De Bremaecker & Swenson 1990;Shen & Stephansson 1995) or a minimum strain-energy density criterion (Sih 1974). As was pointed out by De Bremaecker & Wei (1994) and Wei & De Bremaecker (1995b), these other criteria are incapable of accounting for friction and are therefore not suitable for the problems of closed or mixed-mode fractures. Ignoring this fact may have reinforced the generally held belief that shear cracks cannot propagate in their own plane; this statement is uncertain or wrong (De Bremaecker & Wei 1994).
The ¢rst step in the numerical realization is to calculate the strain energy per unit length perpendicular to the fracture and corresponding to the given fracture with length l (e.g. Hahn 1976), It turns out that *P is important in simulating the upward propagation of £uid-¢lled fractures; it is therefore included in W. The fracture is then enlarged with a de¢ned incremental length dl at one side in a direction r, and the corresponding strain energy is calculated. Usually, dl/l¦0X1. The length-normalized change in strain energy for incremental growth in direction r, the so-called Gri¤th force, is Gr(r){ [W (lzdl){W (l)]/2dl. If possible the fracture will grow in the direction for which Gr is largest. To estimate this direction we ¢rst test three trial directions, one in the planar continuation of the last increment of the fracture and two deviating about a ¢xed angle from that direction (e.g. 10 0 , see Fig. 1). The direction of maximal Gr is then calculated from inverse parabolic interpolation (e.g. Press et al. 1992), leading to a best angle of incremental growth between {10 0 and 10 0 from the direction of the last increment. This scheme works su¤ciently well for our purposes. Whether the direction with maximal Gr is ¢nally accepted and the fracture is enlarged depends on whether Gr exceeds the necessary value to break the rock, which can be estimated by K 2 c (1{l)/(2G) for fractures in homogeneous media and under homogeneous stress, where K c is the fracture toughness, G is the shear modulus and l is the Poisson's ratio. In the following, we often assume a vanishing fracture toughness so that the fracture will always grow when Gr b 0.
It may be observed that the opening dislocation *u N at the fracture tip is zero. Then, the fracture closes and its length is decreased by dl. So, for every iteration or pseudo-time step, a given fracture will increase by dl in the direction of maximal Gri¤th force, remain unchanged or decrease by dl. The two ends of the fracture are treated alternately.
We have included the shear dislocation and shear traction in the formula for W. For simple tensile fractures, and for most £uid-¢lled fracture problems, B S j and *u S j will both be zero or small and will not contribute to W. However, when dealing with £uid-¢lled fractures under deviatoric stress and stress gradients the fractures may propagate in a direction where B S j and *u S j are non-zero. This is explicitly shown below. Under large external shear stresses the shear deformation can become increasingly important in W , and we found some examples for which the £uid-¢lled fracture then starts to propagate as a kind of mixed-mode fracture, in the direction of maximal shear stress. Parts of the opposing fracture surfaces then touch each other. Although this is very interesting, and may have relevance for some geophysical problems, the method is not yet complete enough to consider accurately all aspects of mixedmode fracture propagation. This is because we do not account for sliding friction at the closed parts of the fracture that are moving under shear mode. We have also not included two independent fracture toughnesses for mode I and mode II Figure 1. Sketch showing three trial directions for the incremental growth by dl of the fracture with length l. The ¢nal growth is in the direction of the maximal strain-energy release and is determined from the strain-energy release in the three trial directions and inverse parabolic interpolation. ß 2000 RAS, GJI 141, 623^638 Fluid-¢lled fractures in the Earth fracture propagation (see also Wei & De Bremaecker 1995a;Shen & Stephansson 1995).
The estimation of strain energy needs some additional comments. We will present examples for which ¢ctitious dislocation segments have been used to simulate topographic load on the surface (see also Dahm 1996 for further explanation and veri¢cation). When calculating the strain energy the work done at the free surface segments has to be considered (sum over surface displacements multiplied by loading stress).

PR INCIPLES OF FLUID -FILLED FRACTURE PROPAGATION
The quasi-static growth and propagation of £uid-¢lled fractures is controlled by internal £uid pressures, external elastic stresses and fracture toughness. It is more complicated than for empty tensile cracks since the £uid is balancing its internal pressure distribution immediately after a change in position or in volume occurs. Fracture propagation from one place to another does not usually occur for empty fractures in media under compression, but is possible for £uid-¢lled fractures. While empty fractures usually grow in the direction of maximal compressive stress, £uid-¢lled fractures may grow and propagate in other directions. This section discusses the controlling factors of £uid-¢lled fracture propagation. Because we have magma-¢lled dykes in mind, the parameters have primarily been chosen to demonstrate major e¡ects, and are not necessarily realistic. Applications will be presented later in this paper. Compressive stress is positive-valued in the following.

E¡ect of deviatoric homogeneous stress ¢elds
In the ¢rst examples we concentrate on £uid-¢lled fractures without apparent buoyancy forces and give a comparison to empty tensile fractures. Apparent buoyancy forces will be absent when the £uid has the same density as the surrounding material, or for horizontal fracture propagation (a brief explanation of apparent buoyancy is given in Section 4.1). Fig. 2 shows a comparison of empty and £uid-¢lled fractures of length 2a~5 km under di¡erent stress ¢elds. The rock shear modulus and Poisson's ratio here and below are G~30 GPa and l~0X25, respectively. The £uid has a bulk modulus of K~20 GPa and an initial volume per unit length perpendicular to the cross-section of A 0~1 0 300 m 2 . We assume an initial overpressure for both fractures, under vanishing external stress, of P 0~2 MPa. This leads to a maximal opening in the middle of both fractures of 2X5 m. The opening has an elliptical form identical to the theoretical one; the detailed opening and form of the fractures are not shown in the ¢gure. Under vanishing external stress both fractures are stable, i.e. not growing, since the stress intensity factor K I at their tips is equal to the assumed fracture toughness K c of the rock, K c~KI~P0 na p1 77 MPa m 1a2 . Applying external compression, tension or deviatoric stress alters the shape of the fractures and may lead to a growth of the fractures. This is discussed below.
First, introducing an external isotropic compression of p 1~p3~2 MPa leads to a di¡erent behaviours of the empty and £uid-¢lled fractures (Fig. 2, I). Both fractures are still not growing, but while the empty fracture has completely closed itself, the £uid-¢lled fracture is practically unchanged in shape with a maximal opening of 2X5 m. The compressed £uid within the fracture can sustain the external compression with only a small change in volume. For the empty fracture *u~0 since the external compressive stress exactly compensates the internal overpressure within the fracture. The non-penetrating constraints (inequality constraints) are e¡ective when solving the boundary value problem for the empty fracture.
Applying an external isotropic tension of p 1~p3~{ 2 MPa leads to a two-sided growth of the empty fracture in its own plane (Fig. 2a, II). Due to the tensional stress the fracture opening is larger than for the fracture under zero external stress, leading to an increased stress intensity factor at the tips that exceeds the fracture toughness, K I b K c . Therefore, growth starts. The growth direction is controlled by the orientation of the initial fracture, and the fracture would grow to in¢nity with an increasing maximal opening with increasing length. We stopped calculations after 50 iterations. The £uid-¢lled fracture, in contrast, only propagates a very short distance (less than 0.05 of the fracture length) and is practically stable. Although external tensional stress is applied, the form of the fracture and thus the stress intensity factor at the tips is practically unchanged, due to the large bulk modulus of the £uid. To demonstrate how the £uid-¢lled fracture may grow under more favourable conditions we set the fracture toughness to zero. The fracture then starts to grow in its own plane but stops after a certain length (about 3 km). In Fig. 2(b, II) we have plotted the ¢nal fracture as a dashed line. The pressure in the £uid and the width of the fracture steadily decrease with growing fracture length. The growth stops after a length is reached at which further growth would lead to underpressure within the fracture and thus to closure of the fracture. Note that the situation described above would not occur in reality, since the phase of the £uid would change before a zero or negative pressure was reached. However, the stopping of £uid-¢lled fracture growth due to decreasing internal pressure can occur under compressional conditions, and the simulation above is therefore instructive.
Applying remote stresses of p 1~pyy~2 MPa, p 3~pxx{ 2 MPa demonstrates the e¡ect of deviatoric stress on tensile fracture growth (Fig. 2, III). The fracture faces are not in contact for either the empty or the £uid-¢lled fracture. The fractures have a tendency to grow in the direction of the local maximal compressive stress and to open in the direction of least compressive stress (here the tensile stress). However, recall that we use an energy criterion to estimate the propagation direction. Since the fractures are unfavourably oriented with respect to tensile fracture growth, their growth path is initially more complex. The initial fractures experience shear forces and thus shear dislocations, which leads to a complicated stress ¢eld in the neighbourhood of their tips. The curved growth paths in Fig. 2 conserve approximately the orientations of the maximal compressive stress at the fracture tip during growth. The paths di¡er slightly for the empty and £uid-¢lled fractures since the maximal opening, and thus the tip stress, are di¡erent in both cases. The empty fracture again propagates to in¢nity, while the £uid-¢lled fracture is practically stable for K c~1 77 MPa m 1a2 and grows a ¢nite distance for a vanishing fracture toughness (larger than the plot range in this case).
These simple examples demonstrate that under homogeneous stress, fractures grow at both ends, but do not undergo wholesale displacement from one place to another. The growth length of £uid-¢lled fractures under external load is small since they practically have a constant opening and a changing internal pressure with increasing length.

E¡ect of stress gradients
As is known from theoretical and observational studies of solidi¢ed magma-¢lled dykes (Pollard & Mueller 1976), stress gradients change the form of tensile fractures. As a consequence, such fractures may start to propagate in the direction of stress decrease, even under purely compressional conditions as in the lithosphere. The e¡ect of stress gradients on £uid-¢lled fracture propagation can be demonstrated by the numerical method. While the internal £uid pressure is constant (for example, for horizontal fractures or lateral fracture propagation at the level of neutral buoyancy), the normal stress in the rock and away from the fracture decreases or increases. The balancing of the elastic stresses with the internal constant pressure leads to a tear-drop-shaped fracture. Similar to the theory ¢rst presented by Weertman (1971a), fractures exceeding a critical length begin to close at one end and simultaneously to break the rock at the other end, and in this way may propagate large distances. The real situation can be quite complicated when deviatoric stresses are present in addition to a gradient, or when di¡erent gradients act on the three eigenvalues of the stress tensor. Fig. 3 shows the simulation of £uid-¢lled fractures within a compressional stress ¢eld for which all components increase in the x-direction with the same gradient. In contrast to Fig. 2, both fracture tips move in the same direction and the fracture has an approximately constant length during propagation. The fractures propagate in a direction for which the normal stress on the fracture decreases. Under isotropic compressional stress the propagation direction slowly changes to the negative direction of the stress gradient ({x-direction, Fig. 3 left). The curvature of the propagation path decreases with increasing fracture length (increasing cross-sectional area A 0 ). Fractures with a larger A 0 generate a larger stress disturbance in front of their tips and therefore need a longer propagation length until they propagate into the direction of maximal stress decrease in Fig. 3. It is observed that in this case, as in many of the following examples, the fracture undergoes both a growth in length and a wholesale displacement, resulting from growth at one end and closure at the other.
When adding 10 MPa to p yy and leaving the other stress components unchanged, a deviatoric stress is present with the maximal compressive stress p 1 pointing in the y-direction. Two competitive mechanisms become e¡ective: (1) the fracture has a tendency to grow in the p 1 -or y-direction; and (2) the driving force for fracture propagation is strongest in the direction of maximal stress decrease, i.e. the {x-direction. As a compromise, the fracture path is inclined at about 13 0 against the p 1 -direction in Fig. 3 (middle). This angle increases with decreasing magnitude of p 1 . The fracture length in this case is relatively large since its average opening is relatively small due to a relatively small e¡ective stress gradient along this path.
When adding 10 MPa to p xx instead of p yy , the p 1 -direction points into the direction of maximal stress decrease. Therefore, both mechanisms favour propagation in the {x-direction, and the fracture path in Fig. 3 (right) immediately turns to this direction.
These examples demonstrate that £uid-¢lled fracture propagation can occur under purely compressional stress when ß 2000 RAS, GJI 141, 623^638 stress gradients are present. The driving force is the gradient of the normal stress component along the fracture; that is, fractures cannot propagate along isolines of the normal stress component. The propagation needs neither to be in the p 1 -direction nor to be in the direction of maximal stress decrease, but may be in between the two.

Stress gradients and large deviatoric stresses
The £uid-¢lled fracture propagation can be even more complicated when stress gradients and large deviatoric stresses are present. Then, a type of mixed-mode fracture propagation may occur. This can be demonstrated when using a di¡erent gradient for p xx (30 MPa km {1 ) and p yy (50 MPa km {1 ) in Fig. 3(c) so that the deviatoric stress increases linearly in the {x-direction, while the stress magnitudes decrease in that direction (see Fig. 4). p 1 always points in x-direction. The fracture ¢rst changes its initial direction and propagates into the direction of maximal stress decrease and of p 1 (both in the {x-direction), similar to Fig. 3(c). At x&{7 km the deviatoric stress is large enough to promote propagation into the direction of maximal shear stress, i.e. inclined at 45 0 to the x-axis. Mixed-mode £uid-¢lled fracture propagation may be important for some Earth problems, for example, £uid migration in the stress ¢eld of a subducting slab and possible initiation of deep earthquakes. However, as already mentioned, we think that our approach is not yet ready to simulate mixedmode fracture propagation accurately. Sliding friction, di¡erent material strengths for shear and tensile fracturing, and a possible fracture segmentation and the resulting pressure drop within the £uid are not considered. In the following examples we will avoid mixed-mode fracture propagation for which parts of opposing fracture surfaces are in contact.

E¡ect of density di¡erence between £uid and rock
So far we have assumed that the £uid in the fracture has the same density as the surrounding rock, or that fractures propagate horizontally. When dropping this assumption an apparent stress gradient along the vertical extent of the fracture can occur; this is often called a buoyancy e¡ect. For rocks denser than the £uid, the lithostatic stress decrease in the rock in the direction towards the surface is larger than the hydrostatic stress decrease in the £uid, resulting in an increased overpressure at the top of a vertical fracture. This favours fracture propagation upwards, very similar to the fracture propagation into the direction of maximal external stress decrease. One may speak of apparent Archimedean buoyancy forces acting on the fracture (Weertman 1971a). The buoyancy e¡ect is important for £uid-¢lled fractures in the lithosphere. It can also be modelled with the numerical technique. In the examples above we prescribed a constant £uid pressure on the walls of the fracture (P 0~2 MPa). Now we prescribe a £uid pressure increasing with depth with a gradient of og (o is the £uid density and g is the gravity). The average £uid pressure used is not so important since it is automatically balanced when solving the dislocation problem. Then, every time the fracture growth has a non-zero vertical component *z, we consider the change of the internal pressure at the new fracture tip. The pressure at the new fracture tip is estimated by adding the Figure 3. Simulation of £uid-¢lled fracture propagation in compressional stress ¢elds for which all components increase at a rate of 50 MPa km {1 in the x-direction. On the left p 1~p3~8 00 MPa at x~0. In the middle 10 MPa has been added to p yy , and on the right to p xx . The initial and ¢nal fractures after 100 iterations are indicated as bold lines; the initial fracture is near the origin (lower right corner). K c~0 has been assumed, and the initial fracture has not reached its equilibrium length (that is, the stress intensity factors at both tips are initially overcritical). The fracture propagation path is given as a dashed line. A 0 is the initial cross-sectional area of the fracture; in the middle and on the right A 0~2 000 m 2 was chosen. pressure value from the former tip and *P~og*z. In addition to the pressure gradient in the £uid, a vertical gradient in the external stress is considered, for example, resulting from the gravitational stress in the lithosphere. Fig. 5 shows simulations of buoyancy-driven fractures from the lithospheric mantle into the lower and upper crust, that is, in a two-layer over half-space model. We used rock densities of 3300, 2900 and 2600 kg m {3 for the mantle, lower crust and upper crust, respectively (see Fig. 5). The £uid density is 2950 kg m {3 (e.g. dense tholeiitic magma). The apparent buoyancy force is directed upwards in the mantle and downwards in the crust; the crust^mantle boundary is the so-called level of neutral buoyancy (LNB). The initial fracture, or initial magma reservoir, is at 100 km depth. After initiation, the fracture propagates upwards by breaking the rock at its top and simultaneously closing itself at its tail. For a crust under zero deviatoric stress or under horizontal extension, the fracture`overshoots' the LNB without a change in the vertical propagation direction. The propagation stops somewhere in the upper crust, depending on A 0 and on the speci¢ed density model (Fig. 5a). All dykes still have their bottom end in the upper mantle after propagation has stopped.
For a crust under horizontal compression the propagation path and the ¢nal stable position of the fracture depend strongly on A 0 (Fig. 5b). Along the mantle path, for A 0~5 |10 5 m 2 , the length and the maximal opening of the fracture are 30 km and about 23 m, respectively. The maximal overpressure at the upper tip is about 78 MPa. For A 0~5 |10 3 m 2 the corresponding values are 6 km, about 1 m and about 15 MPa. Using the large initial volume of A 0~5 |10 5 m 2 the fracture is able to overcome the compressional stress barrier in the crust and to reach the free surface. For smaller values of A 0 the fracture propagates in the horizontal direction and stops when the driving force has fallen below a certain limit. Most of the length of each fracture is then horizontal (Fig. 5b). The observed small downward propagation for one subhorizontal fracture in Fig. 5(b) is most probably an e¡ect of the stress ¢eld generated by the fracture itself. The sill-like bodies that form are always above the LNB, and their depth is larger for fractures with a smaller A 0 .
Magmatic underplating has been frequently invoked as a mechanism of crustal growth, where sill-like intrusions in both the upper mantle and the lower crust are discussed (e.g. Behrendt et al. 1990;Deemer & Hurich 1994). Often it is suggested that underplating occurs in extensional regimes. Our examples suggest that dense ma¢c intrusions are most probably arrested in the lower crust, and that sheeted dyke and sill complexes may be generated there under extensional and compressional tectonic stresses, respectively. For instance, the highly re£ective lower crust typically observed in Hercynian and Caledonian consolidated continental crust (e.g. Matthews & Cheadle 1986;Lu« schen et al. 1987) represents layered ma¢c intrusions that may have been formed by dyke intrusions from the mantle and sill formation during the compressional orogenesis. This agrees with the view of Wenzel & Sandmeier (1992), who suggested a`dry' and layered basaltic intrusive structure of the highly re£ective lower crust beneath the Black Forest, SW Germany, to explain joint P-and S-wave observations.
Another noteworthy result of our numerical simulations, not documented in Fig. 5, is that horizontal compression is able to hinder the rise of £uid-¢lled fractures even under positive apparent buoyancy forces, and may cause the formation of sills containing a £uid less dense than the surrounding rock. Thus, repeated intrusions can form potentially unstable £uid reservoirs in compressional tectonic regimes. The reservoir may, after reaching a critical volume, be able to fracture the rock again and overcome the stress barrier in a single event.
Note that the results apply to 2-D problems (plane strain); the 3-D situation may be more complex.

E¡ect of topographic load on fracture propagation
Surface topography is quite common in regions where enhanced £uid ascent and accumulation is observed, for example, beneath volcanic rifts, volcanoes and back-arc volcanic chains, as well as £uid accumulation beneath orogenic belts and continental graben systems. Magma and aqueous £uids are a¡ected. Surface topography has a far-reaching impact on the stress ¢eld and thus in£uences the propagation paths of £uid-¢lled fractures and the regions of preferred £uid accumulation.
For instance, the strong in£uence of gravitational stress on the propagation of dykes at the Hawaiian volcanic rifts has already been noted by Fiske & Jackson (1972). Due to gravitational stress, dykes propagate laterally up to 120 km along the axis of pre-existing rifts from the volcanic centres of the rift to its tips. Fiske & Jackson (1972) also demonstrated this behaviour with laboratory experiments using solidi¢ed gelatin rifts' and water injections.
A simple method to study the topographic e¡ect on fracture propagation is to apply a topographic load at the surface of a homogeneous elastic half-space. For a non-interactive problem the topographic load may be calculated from the gravity stress generated by the overlying rocks. In this model the elasticity of the substratum is unimportant to the overall response of the half-space. These are situations where the substratum is either very £exible or very rigid when compared with the half-space elasticity. Interactive problems consider the coupling between the overlying rock and the half-space (e.g. Davis & Selvadurai 1996), which is much more complicated and not treated here. Bending of the lithospheric plate may be an additional unconsidered aspect of large-scale problems.
Before simulating fracture propagation we discuss analytical stress ¢eld solutions of a uniform and a pyramid-shaped strip load of length 2b acting on an elastic half-space. The solutions are based on Flamant's fundamental solution for a line load on an elastic half-space and are brie£y described in Appendix A.  uniform and pyramid-shaped strip loads. Both solutions depend only on the problem geometry and the maximal loading stress P 0 , and are independent of elastic constants. At x~z~0 the maximal (p 1 ) and minimal stresses (p 3 ) are both compressive and equal to P 0 . p 1 decreases to a quarter of the maximal loading stress P 0 at depths of about 5 and 2.5 half-widths b of the topographic load for the uniform and pyramid-shaped loads, respectively. The decrease of p 3 is larger in both cases. The stress perturbation at a given distance is larger in the vertical than in the horizontal direction. The gradients of p 1 and p 3 as seen in Fig. 6(b) would promote downwards fracture propagation diverging from the centre of the topographic load. The observed inclination of the p 1 -axes promotes fracture propagation towards or away from the centre of the topographic load. As is shown below, the apparent buoyancy forces acting on fractures containing a £uid less dense than the surrounding rock lead to an upward propagation of the fractures that is controlled by both the buoyancy forces and the direction of the p 1 -axes. We show that the increase of p 1 and p 3 in Fig. 6(b) towards the centre of the loading at the surface is able to prevent £uid-¢lled fractures from actually reaching the surface in some cases.
Another interesting aspect is when in addition to topographic loading, tectonic stress is applied (Figs 6c and d). An extensional stress component weakens the in£uence of topographic load on fracture propagation; that is, the fracture propagation paths will deviate less from the vertical upwards direction. In contrast, under compressive tectonic regimes the surface topography is able to rotate the formerly horizontal p 1 -axis to a vertical direction beneath the strip load, and may thus generate a stress channel that may enable £uid-¢lled fractures to propagate upwards and reach the surface.
The e¡ect of topographic load on fracture propagation is demonstrated in Figs 7^10. The free surface and the topographic load were simulated using 82 segments, respectively stress-free and representing appropriate normal stresses in the range {9¦x/b¦9. Such a large horizontal extension of the model beyond the plotted range in Figs 7^10 was found to be necessary by cross-checking the numerical fracturefree solution against the analytical solution. In addition, a lithostatic isotropic stress linearly increasing with depth at a rate of 7.5 P 0 /b has been superposed. This stress gradient is of minor interest and therefore not further discussed. The initial positions of the fractures were at z/b~{4, and at various horizontal positions in the range {4¦x/b¦4. In addition to buoyancy forces, gradients and principal directions of the background stress, the fracture propagation is in£uenced by the initial volume, the length of the fractures and the overpressure within the fractures resulting from the density di¡erence between host rock and £uid. A typical length can be determined from the analytical expression of a vertical fracture in a homogeneous space, (e.g. eq. 4 of Weertman 1980). The overpressure *P~*ogl gives a rough estimate of the pressure at the tip of the vertical fracture, and is closely related to the upward fracture-driving force. We study the fracture propagation as a function of *P/P 0 and l/b, where P 0 and b are the head pressure and the half-length of the strip load at the surface, respectively. Fig. 7 gives two examples of £uid-¢lled fracture propagation in a homogeneous half-space under lithostatic stress with a On the left (a) the apparent buoyancy force is smaller (*P/P 0~0 X185) than on the right (b, *P/P 0~1 X11); l/b~0X25 in both cases. Interaction between fractures is not considered. ß 2000 RAS, GJI 141, 623^638 pyramid-shaped load applied at the surface. If the theoretical overpressure *P is only about 20 per cent of the head pressure P 0 of the load (Fig. 7a), the fractures at di¡erent horizontal positions at depth follow the direction of the maximal compressive principal stress p 1 and the propagation paths converge to the centre of the load at the surface. All fractures reach a depth z~0. For instance, if b&30 km and P 0 &120 MPa, similar to the topographic load of Kilauea volcano, Hawaii, magma-¢lled dykes with a length of 7.5 km and a *o of about 300 kg m {3 are predicted to accumulate beneath the summit of the volcano. Dykes at a horizontal distance and depth 4b are still a¡ected by the topographic load.
The deep structure beneath mid-ocean ridges, for example, as found by the mantle electromagnetic and tomographic experiments (MELT) at the East Paci¢c Rise (e.g. Forsyth et al. 1998), may be another application. A broad asymmetric pyramid-shaped region of low seismic velocities extending to a depth of 150 km and ranging from 400 km west to 250 km east of the ridge is interpreted to be the primary melt production region beneath the ridge. The question is which is the mechanism guiding the melt from this broad region to the narrow ridge axis? Our work suggests that magma is ascending via propagating fractures and that the gravitational stresses of the ridge guide the melt to the ridge axis. Fig. 7(b) shows that the propagation paths of £uid-¢lled fractures are not only controlled by the direction of p 1 but also by the apparent buoyancy force *P/P 0 . Using a *P/P 0 of 1.1 leads to a deviation of propagation paths from the local p 1 direction for dykes far from the centre of the loading, where magnitudes of the background stress are smaller. The typical distance between neighbouring volcanoes may be controlled by the competition between the inhomogeneous stress ¢eld produced by their own weight and the buoyancy of the fractures. A somehow similar model has been proposed by Hieronymus & Bercovici (1999) to explain the alternating hotspot islands in the Hawaiian chain. However, this model does not take into account the large e¡ect of the inhomogeneous stress ¢eld on the direction of fracture propagation.
Adding a compressive tectonic stress of p xx~P0 /12 alters the propagation paths in such a way that the £uid-¢lled fractures use only a small channel under the strip load to propagate upwards (Fig. 8). A typical situation where topographic load is Figure 8. The same model as in Fig. 7 but with an additional tectonic compressive stress of p xx~P0 /12. The fractures are characterized by *P/P 0~0 X185, and by l/b~0X25 (a) and l/b~0X39 (b). Note that an increase of l/b and a constant *P/P 0 are accompanied by a decrease in the density di¡erence between rock and £uid. Figure 9. Comparison of the opening dislocation (discretized fracture opening) of a fracture at x~0 and an upper tip depth at z/b~{0X25 with theoretical fracture lengths of l/b~0X25 (dashed line) and l/b~0X39 (solid line) and *P/P 0~0 X185. For b~30 km, P 0~1 20 MPa, G~30 GPa and l~0X25, the opening dislocation is given in metres. The discretization of fractures used can be seen from the segment length with constant opening. ß 2000 RAS, GJI 141, 623^638

632
T. Dahm acting in addition to compressive horizontal stress is beneath an island arc at a destructive plate margin or beneath an active continental margin such as the Andean Cordillera of South America. The observed topography at 25 0 S at the eastern continental margin of the South American Plate is up to 12 km when measured from the trench and has a half-width of about 250 km. The regional stress ¢eld in continental South America is dominated by E^W compression (Zoback 1992). Such a major topographic feature should signi¢cantly a¡ect the propagation of £uid-¢lled fractures even at a depth of 200 km, where tomographic experiments and seismicity indicate that large amounts of £uids leave the subducting slab and migrate into the overlying mantle wedge (e.g. Haberland & Rietbrock 1999). However, due to the viscous deformation of the mantle wedge and the subduction of the Nazca Plate, the stress ¢eld at 200 km depth may be more complex than that estimated from the topographic load alone; this is further discussed in the next section. It is interesting that dykes with a typical length of l/b~0X25 and *P/P 0~0 X185 reach the surface, while those with l/b~0X39 and *P/P 0~0 X185 do not (Fig. 8); the latter stop propagating upwards at a depth of about b/4. As mentioned above, the decrease in the gradient of p 3 towards the surface beneath the strip load may prevent fractures from reaching the surface, although the £uid density is less than the rock density. This di¡erence in fracture propagation would not be easy to predict from simple analytical full-space models, for example, from looking only at the level of neutral buoyancy. Fig. 9 compares the opening dislocations of both types of fractures with a theoretical l/b~0X39 and 0X25 at x~0 and an upper tip depth at z/b~{0X25. The longer fracture is stable in this position while the smaller one further propagates upwards to the surface. The longer fracture in Fig. 9 has an approximately tear-dropshaped opening with the narrow end at the top, indicating that the apparent buoyancy force is nearly zero in this position since fracture propagation stops when further growth is accompanied by an increase of strain energy or when the fracture does not open further as its length grows.
Increasing the tectonic compressive stress to P 0 /6 leads to an interesting result. Fluid-¢lled fractures with *P/P 0~0 X185 and l/b~0X25 do not propagate upwards; instead, they accumulate in horizontal sills beneath the centre of the strip load at a depth range between {4 and {2X5 in our examples (Fig. 10a).
When the theoretical overpressure *P/P 0 is 1.1, the snapshot after 220 pseudo-time steps in Fig. 10(b) shows a complex situation. Fractures initiated at a depth of z/b~{4, outside the horizontal range of the strip load, are able to propagate directly upwards when entering the stress channel that has a non-horizontal p 1 -direction. Fractures initiated at jx/bj`1 and z/b~{4 ¢rst form a horizontal sill in the depth region where the stress component from the strip load is relatively small. Later they stop propagating in one direction and start to propagate in the opposite direction but now upwards to ¢nally reach the stress channel where p 1 is steeply inclined and where the stresses from the strip load are relatively large. It is clear that such a complex and long propagation path may be altered or disturbed in reality, for example, when magma-¢lled fractures solidify at depth or when propagation paths cross pre-existing fractures or zones of weakness.
It should be noted that this section has dealt with a single fracture propagation in a half-space. More complex cases could be similarly investigated.

Fluid ascent in the mantle wedge above a subducting slab
The motion of a large plane slab that sinks beneath an island arcs controls the £ow of and thus the stress within the mantle wedge if other forces are absent. This £ow pattern is in a ¢rstorder approximation similar to the stationary corner £ow that has been analytically derived (e.g. Batchelor 1967;McKenzie 1969) assuming a constant viscosity g in the mantle wedge, a slab dipping with an angle # S , and a subducting velocity o. The shear stress p r# is given by (compressive stress positive) Figure 10. The same model as in Fig. 7 but with an additional tectonic compressive stress of p xx~P0 /6. The fractures are characterized by l/b~0X25, and by *P/P 0~0 X185 (a) and *P/P 0~1 X1 (b). The simulation was stopped after 220 iterations, and positions of fractures represent a corresponding snapshot.
where the origin of the cylindrical coordinate system is at the corner of the mantle wedge, and r and # give the radial and tangential directions, respectively. Note that # is negativevalued since the z-axis is pointing upwards. p rr and p ## are both zero. The largest shear stresses are thus parallel to the r~const. and h~const. directions, and stress components are proportional to og and decrease with reciprocal distance r {1 . At the line the shear stress is zero. The maximal and minimal principal deviatoric stresses are compressive and tensional, respectively, with their strongest values at the surface of the slab and the lithosphere and zero values at #~# 0 . Thus, the decrease of the least compressive stress p 3 (isotropic lithostatic plus deviatoric stress) along a line r~const. in the upward direction is smaller than the lithostatic stress for # S`#`#0 and larger for # 0`#`0 0 . A £uid-¢lled fracture may experience a positive apparent buoyancy force outside and in the upper part of the wedge, but a negative buoyancy force in the lower part of the wedge. Eq. (5) is assumed to be a good approximation to the stress ¢eld in the mantle wedge for r §50 km (McKenzie 1969).
It is generally agreed that the process of magma generation in a subduction zone is a multistage multisource process (e.g. Wilson 1989). Slab-derived aqueous £uids and partial melts propagate in the mantle wedge above the slab and lower the solidus there, thereby promoting partial melting. They may be regarded as a catalyst for the bulk of subduction zone magmatism rather than a primary source. The transport mechanism of the slab-derived £uids may be fracture propagation, and will then be controlled by apparent buoyancy forces and the stress gradients and principal axis directions in the inhomogeneous stress ¢eld of the mantle wedge. Fluid ascent via fracturing within the mantle wedge above the slab has been proposed in several papers (e.g. Davies 1999;Davies & Stevenson 1992), but has never been quantitatively tested.
In Fig. 11 we simulate the magma-¢lled dyke propagation in the inhomogeneous stress ¢eld of the mantle wedge above a subducting slab (dip angle # S is {45 0 , g is 0X5|10 20 Pa s and o is 8 and 0.8 cm yr {1 in Figs 11a and b, respectively). A lithospheric thickness of 50 km and a rock density of 3300 kg m {3 are assumed. The magma has a density of 2950 kg m {3 and the dykes have a vertical length of 14 km in the shear-stress-free region below and above the mantle wedge corresponding to A 0~5 |10 4 m 2 . The propagation increment is 1 km.
In our model the shear stress is non-zero only within the mantle wedge. At the slab surface in Fig. 11 (#&# S~{ 45 0 ) the compressional principal stress p 1 is horizontally oriented. Fractures propagating with a small apparent buoyancy force and with a short length tend to form sills at this positions. Further inside the wedge, p 1 rotates in such a way that the branch pointing upwards is directed away from the slab. When the line #~# 0 is passed, where the shear stress is zero, the upward-pointing branch of p 1 is directed towards smaller distances r so that fracture propagation is promoted upwards and towards the corner of the wedge. In addition, factors in£uencing fracture propagation are the apparent buoyancy forces, the stress gradients and the self-disturbed stress ¢eld in front of the fracture tip. The in£uence of all of these components can be recognized in Figs 11^13 in di¡erent regions in the wedge.
Apparent buoyancy forces depend on *o and on the vertical length of the fractures and are strongest in the z-direction (upwards). Their in£uence is large if the deviatoric stress is small, which leads to nearly vertical propagation paths, for example, far from the corner of the mantle wedge or for small values of og (Fig. 11b). At small distances r the orientation of the principal stress is the major in£uence on fracture propagation paths. Therefore, the fracture starting at the slab at a horizontal distance of only 30 km from the corner nearly follows the local upward-pointing branch of p 1 through the wedge and reaches the lithosphere at a horizontal distance of 70 km (Fig. 11a). The stress ¢eld of the mantle wedge leads to a horizontal o¡set *x away from the wedge corner for magma dykes ascending from the surface of the slab (Fig. 11a). As is shown below, the same may occur for water ascending from the surface of the slab.
How does the situation change for a slab subducting at a steeper or more shallow angle than 45 0 ? For instance, a large dip angle (about 60 0 ) is observed at 100 km depth for the slabs at the Lesser Antilles, Marianas and New Hebrides. Shallowdipping subduction of 20 0 and less can be found for the Nazca plate under South America. The absence of active volcanism beneath anomalously shallow-dipping segments (`10 0 ) of the Nazca plate is striking and has been attributed to the displacement of the asthenospheric mantle and the direct superposition of the two lithospheric plates (Barazangi & Isacks 1979). In Fig. 12 subduction angles # S of {20 0 and {60 0 are tested. For # S~{ 20 0 the deviatoric stress is a bit larger than for # S~{ 45 0 if og is unchanged; for # S~{ 60 0 it is smaller. A second change compared to the # S~{ 45 0 case is that the upward-pointing branch of p 1 at the slab surface is directed upwards and away from or towards the slab for # S~{ 20 0 and {60 0 , respectively. Thus, for # S~{ 60 0 (in general for j# S j b 45 0 ) there is a possibility that £uid-¢lled fractures propagate upwards and towards or along the slab surface and do not traverse a large part of the mantle wedge, as is observed in the simulation of magma-¢lled dyke propagation in Fig. 12(b), except for fractures at large horizontal distances from the wedge corner ( §140 km), where the deviatoric stresses are relatively small. (In fact, these fractures started at z~{200 km within the wedge.) For the shallow-dipping slab (Fig. 12a) the propagation paths are similar to the # S~{ 45 0 case for distances x b 90 km. However, the horizontal o¡set experienced by a fracture that traverses the mantle wedge is slightly smaller than for # S~{ 45 0 , mainly because the wedge is narrower. Note in particular that dykes that enter the mantle wedge at distances x¦80 km do not reach the lithospheric layer but form horizontal reservoirs (sills) 3^5 km below the lithosphere (Fig. 12a). This is due to the strong curvature of the p 1 direction and the relatively large deviatoric stresses in the narrow wedge in this region. The dykes do not follow the sharp upward bend of the p 1 axis there and prefer to propagate a small distance downwards in the p 1 direction, until propagation stops due to the lack of buoyancy (i.e. most of the fracture is horizontal). The dykes are relatively long in the inhomogeneous stress ¢eld of the mantle and therefore have an apparent inertia due to the self-made perturbed stress ¢eld at their tips. This prevents them from making sharp bends. Their length is increased due to the smaller apparent buoyancy force they experience when Figure 12. Simulations similar to Fig. 11(a), but with slab dip angles of {20 0 (a) and {60 0 (b). Figure 13. Simulation of water-¢lled fracture propagation. The geometry and boundary conditions are the same as in Fig. 11. The water-¢lled fractures have a smaller net apparent buoyancy force than the magma-¢lled fractures in Fig. 11, mainly due to their short length in this simulation. The fracture length is in the range of 400 m and can barely be resolved at the ends of the propagation paths. ß 2000 RAS, GJI 141, 623^638 they propagate in a nearly horizontal direction. This result gives rise to the question whether a narrow mantle wedge and a fast convergence velocity may prevent magma-¢lled dykes from traversing the wedge, and therefore lead to a cessation of volcanism as observed at the shallow-dipping segments of the Nazca plate.
Another question is whether the propagation path in the mantle wedge and the horizontal o¡set have an in£uence on the position of volcanic chains and volcanic islands above subducting slabs. A hint may come from attenuation tomographic studies beneath volcanic arcs. Haberland & Rietbrock (1999) analysed seismic waveforms recorded at a local seismic network in the Central Andes. Their tomographic images indicate that the magmatic source region beneath the Central Andes extends down to the surface of the slab (between 100 and greater than 170 km). The deep source region is more than 100 km wide. The regions of high attenuation indicate that the £uids migrate from the broad deep source region near the slab to crustal depths and seem to be stored in a volume beneath the volcanic chain. These results ¢t su¤ciently well with the propagation paths predicted in this paper.
As mentioned before, a signi¢cant amount of slab-derived aqueous £uid may propagate into the overlying mantle wedge. Water-¢lled fractures di¡er from magma-¢lled ones because the apparent buoyancy force per unit vertical length is much larger for water than for magma. We assume a density of 1000 kg m {3 for water. As a consequence, water-¢lled fractures start to propagate when their lengths are much smaller than that necessary for magma-¢lled dykes to propagate. Their lengths are of the order of several tens of metres. We simulated water-¢lled fractures with a length of 400 m when propagating vertically in the shear-stress-free region of the mantle; this corresponds to A 0~2 0 m 2 . The propagation increment length was 10 m. Note that the net apparent buoyancy forces of the magma-¢lled fractures in Figs 11 and 12 are larger than for the £uid-¢lled fractures below due to the relatively short lengths of the water-¢lled dykes. The values have been chosen to be within a plausible range; however, as will be shown below, the results depend critically on these values. This should be kept in mind. Fig. 13 shows a simulation for #~{45 0 , g~0X5|10 20 Pa s and o~8 and 0.8 cm yr {1 , identical to the situation shown in Fig. 11. In the case of the fast-moving slab (Fig. 13a) the water-¢lled fractures often do not propagate a large distance in the mantle wedge, and do not reach the line of zero shear stress in the mantle wedge at which the magnitude of the apparent buoyancy force changes. The cessation of fracture propagation is a result of the reduced gradient of p 3 in the lower part of the mantle wedge compared to the gradient outside and in the upper part of the wedge. The apparent buoyancy force vanishes along a path in the lower part of the wedge in the direction of p 1 , so that the propagation spontaneously stops. The example thus shows that the lower part of the mantle wedge has the potential to stop fracture propagation, although the density of the £uid is smaller than the density of the mantle rock. This feature is not typical of water-¢lled fractures alone; it depends on the magnitude of the apparent buoyancy force in the wedge, and thus on the density di¡erence between £uid and mantle rock, on the length of the fractures, on the slab subduction velocity and on the viscosity of mantle rock.
If og is reduced by about a factor of 10 ( Fig. 13b), the water-¢lled fractures are able to propagate through the whole mantle wedge and reach the lithospheric layer.
As a last remark we mention that the real situation in the mantle wedge is probably more complex than that presented here, since water and also magma may interact with mantle material, which possibly leads to the cessation of fracture propagation.

CONCLUSIONS
Not only can £uid-¢lled fractures grow at their tips, but they can also propagate from one position to another, because they may fracture the rock at one end, and they simultaneously close at the other end. In this way they may propagate large distances with a nearly constant length and £uid volume (£uid mass). Although this has long be known in geophysics it has not always been taken into account in models of large-scale £uid transport in the Earth.
Driving forces for fracture propagation are stress gradients in the rock and apparent buoyancy forces resulting from density di¡erences between £uid and rock. Contrary to commonly held beliefs, the fracture growth and propagation directions are not controlled only by the direction of the principal stresses, but also by tectonic stress gradients, apparent buoyancy forces and the length and thus the £uid mass of the fractures themselves. It is not easy to predict which of these factors is the major in£uence for a given situation, as the examples presented in the paper make clear.
Some applications to geological problems are listed below.
(1) Fluid-¢lled propagating dykes overshoot the level of neutral buoyancy (LNB), e.g. the Moho, and may form sheeted sills under compressional tectonic stresses in the layer above the LNB. The depths of the sills depend on the stress and the initial volumes of the dykes. Dykes with a very large initial volume are able to break the compressional stress barrier and the layer of negative buoyancy, and reach the free surface.
(2) Topographic load has a far-reaching impact on £uid-¢lled fracture propagation. We found that typically the fractures from depth propagate towards the centre of the topographic load (towards the ridge). Under compressive tectonic stress the topographic load can produce a stress ¢eld channel that enables upward migration of £uid-¢lled fractures in an environment where the fractures would otherwise form horizontal sills. Fracture propagation within the catchment area of a gravitational load may explain how magma from a broad source region at depth beneath volcanoes is guided to the centre of the volcano or the rift, although the sideways or upward propagation path to the surface is often shorter. For instance, shallow crustal magma chambers and intrusive complexes are often formed beneath the highest topography of a volcano.
(3) Our contribution to the understanding of the complex system of melt generation and transport beneath convergent margins is to predict how £uid-¢lled fractures propagate through the inhomogeneous stress ¢eld of the asthenospheric mantle wedge above a subducting slab. We found that propagation paths can be complicated; they depend on the dip angle of the slab, the distance from the wedge corner, the factor og (subduction velocity times mantle viscosity) and the apparent buoyancy force of the £uid-¢lled fractures. For dip angles larger than 45 0 , £uid-¢lled fractures are very likely to propagate upwards parallel to the slab in a small region from the slab surface. For shallow dip angles (e.g. 20 0 ) and near to the wedge corner, our simulated dykes enter the mantle wedge but do not reach the lithospheric layer; instead, £uids are trapped in subhorizontal curved sills. Similarly, active volcanism is absent at the shallow-dipping segments of the Nazca plate in the Andean Cordillera region. When shear stresses in the mantle are relatively large near to the surface of the slab (for example, as a result of a fast moving slab), the £uid-¢lled fractures may stop within the lower half of the wedge if their apparent buoyancy force is not very large, that is, if their length does not signi¢cantly exceed the critical length necessary for buoyancydriven fracture propagation. However, the fractures that propagate from the slab through the mantle wedge and reach the lithosphere have a curved propagation path and usually experience a net horizontal o¡set between 10 and more than 100 km from the wedge corner.
The displacement discontinuity method can explain the wholesale migration of fractures due to their opening at one end and closing at the other. It also shows the importance of inhomogeneous stress ¢elds; such phenomena cannot be taken into account by analytical methods.