Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.geol.msu.ru/deps/petro/Perchuk_publ/gerya_tecton02.pdf Äàòà èçìåíåíèÿ: Mon Apr 22 18:43:11 2013 Äàòà èíäåêñèðîâàíèÿ: Thu Feb 27 23:30:13 2014 Êîäèðîâêà: Windows-1251

TECTONICS, VOL. 21, NO. 6, 1056, doi:10.1029/2002TC001406, 2002

Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation
Taras V. Gerya
1

Institute of Experimental Mineralogy, Russian Academy of Sciences, Chernogolovka, Moscow, Russia

Bernhard Stockhert Å
Institut fur Geologie, Mineralogie und Geophysik, Ruhr-Universitat Bochum, Bochum, Germany Å Å

Alexey L. Perchuk
Institute for Ore Deposits Geology, Petrology, Mineralogy, and Geochemistry, Russian Academy of Sciences, Moscow, Russia

Received 2 May 2001; revised 5 August 2002; accepted 27 August 2002; published XX Month 2002.

[1]

High-p r essure metamo rphic r ock s p r ovide evidence that in subduction zones material can return from depths of more than 100 km to the surface. The pressure-temperature paths recorded by these rocks are variable, mostly revealing cooling during decompression, while the time constraints are generally narrow and indicate that the exhumation rates can be on the order of plate velocities. As such, subduction cannot be considered as a single pass process; instead, return flow of a considerable portion of crustal and upper mantle material must be accounted for. Our numerical simulations provide insight into the self-or g anizing lar ge-scale flow p atterns a nd temperature field of subduction zones, primarily controlled by rheology, phase transformations, fluid budget, and heat transfer, which are all interrelated. They show the development of a subduction channel with forced return flow of low-viscosity material and progressive widening by hydration of the mantle wedge. The large-scale structures and the array of pre ssu re -temp e r a tu re pa th s o bta i n e d b y t he se simulations favorably compare to the record of natural rocks and the structure of high-pressure INDEX TERMS: 3210 Mathematical metamorphic areas.
Ge ophysics: Mode ling; 3660 M i neralogy and Petrology: Metamorphic petrology; 8150 Tectonophysics: Evolution of the Earth: Plate boundary--general (3040); 8045 Structural Geology: Role of fluids; 5120 Physical Properties of Rocks: Plasticity, diffusion, and creep; KEYWORDS: numerical modeling, subduction zones, exhumation of high-pressure rocks, serpentinized mantle, PT paths, rheology. Citation: Gerya, T. V., B. Stockhert, and A. L. Å Perchuk, Exhumation of high-pressure metamorphic rocks in a

subduction channel: A numerical simulation, Tectonics, 21(6), 1056, doi:10.1029/2002TC001406, 2002.

1. Introduction
[2] Subduction zones are the sites of recycling of aging oceanic lithosphere into the mantle. While the general position of the downgoing slab is revealed by Benioff zone seismicity [e.g., Green and Houston, 1995; Ruff and Tichelaar, 1996; Kirby et al., 1996], by seismic tomography [e.g., Spakman, 1991; Karason and van der Hilst, 2001], and by a combination of various advanced geophysical techniques [e.g., Yuan et al., 2000] and the temperature field of subduction zones has been simulated based on various assumptions [e.g., Turcotte and Schubert, 1982; van den Beukel and Wortel, 1988; Peacock, 1990a, 1996] for steady state subduction, the fine-scale structure of the subduction zone is at the limit of or beyond the spatial resolution of geophysical methods. However, an ever increasing number of high-pressure and ultrahigh-pressure metamorphic rocks have been discovered in the last two decades, which according to their unique P/T ratios - must have been exhumed from ancient subduction zones. These rocks prove that material returns to shallow crustal levels, and eventually to the surface, from depths of more than 100 km. This happens either continuously or during certain episodes in the lifetime of a subduction zone, with rates corresponding to plate velocity [e.g., Gebauer et al., 1997; Perchuk et al., 1998; Amato et al., 1999; Rubatto and Hermann, 2001]. The record of these rocks reveals information on the pressuretemperature-time paths followed [e.g., Carswell and Zhang, 1999; Ernst, 1999, 2001; Ernst et al., 1995; Ernst and Liou, 1999] and on the state of stress and the mechanisms of deformation [e.g., Stockhert, 2002] at depths and on length Å scales, which are otherwise not accessible for direct analysis. Numerical simulations, with boundary conditions chosen according to the constraints provided by the record of the natural rocks, can then provide insight into the large-scale flow patterns and into their sensitivity to the prescribed parameters. Our simulations are designed in such a way to

1 Now at Institut fur Geologie, Mineralogie und Geophysik, RuhrÅ Universitat Bochum, Bochum, Germany. Å

Copyright 2002 by the American Geophysical Union. 0278-7407/02/2002TC001406$12.00

6-1


6-2

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

allow for self-organization of the subduction zone, i.e., starting with a simple initial configuration.

2. Simulation of Subduction Zone Processes
[3] The principal features of the subduction process have been previously addressed in numerous analytical studies, numerical simulations, and analogue modeling. Without attempting even a very short review of thermal and kinematic [e.g., Peacock, 1996], fully dynamic [e.g., Schmeling et al., 1999] and combined [e.g., Davies and Stevenson, 1992] subduction related numerical modeling, we refer here instead to some recent examples of numerical experiments providing insight into subduction zone processes of relevance to our study. The relation between the thermal and the petrological structure of subduction zones was examined by, e.g., Peacock [1996], Hacker [1996], Allemand and Lardeaux [1997], and Guinchi and Ricard [1999]. Kincaid and Sacks [1997] found that the thermal structure of the mantle wedge of the overriding plate is controlled by forced asthenospheric flow. Schmeling et al. [1999] investigated the influence of metastable preservation of olivine at the mantle transition zone on the subduction rate. Crustal-scale episodic accretion at convergent plate boundaries was examined by Ellis et al. [1999]; van Hunen et al. [2000] simulated the processes involved in shallow subduction beneath an overriding continental lithosphere. RegenauerLieb et al. [2001] proposed a viable mechanism for the initiation of subduction by sedimentary loading at a continental margin. [4] There is little doubt about the deep burial of highpressure metamorphic rocks being due to subduction with the downgoing slab [e.g., Ernst, 1977; Chopin, 1984; Sobolev and Shatsky, 1990; Schreyer, 1995; Harley and Carswell, 1995]. However, the mechanisms of their exhumation remain a matter of debate and a number of feasible models have been proposed [e.g., Platt, 1993; Hacker and Peacock, 1995; Maruyama et al., 1996; Ring et al., 1999]. The characteristic length and time scales recently identified for the exhumation of high- and particularly ultrahighpressure metamorphic rocks [Gebauer et al., 1997; Amato et al., 1999; Perchuk et al., 1998; Perchuk and Philippot, 2000a; Rubatto and Hermann, 2001] pose further constraints and seem to agree with a corner flow model. In this model, originally proposed by Hsu [1971], and further developed by Cloos [1982], Cloos and Shreve [1988a, 1988b], and Shreve and Cloos [1986], exhumation of high- and ultrahigh-pressure metamorphic crustal slices at rates on the order of plate velocity is driven by forced flow in a wedge-shaped subduction channel. The observation that many exhumed highand ultrahigh-pressure metamorphic rocks have remained undeformed, or are exclusively deformed by some kind of diffusion creep in the presence of a hydrous fluid at very low stress, implying a low bulk viscosity of the material in the subduction channel, is taken to support the concept of forced flow [Stockhert, 2002]. Furthermore, using a one-dimenÅ sional analytical model, Gerya and Stockhert [2002] have Å demonstrated that the exhumation rate in a subduction channel is higher for a Newtonian rheology compared to a

power law rheology, and that the highest rate of exhumation correlates with the lowest shear deformation in the central portion of the channel. [5] In contrast to previous studies, here we investigate the self-organizing evolution of the accretionary wedge and the subduction channel. This implies that the geometry of the accretionary wedge and the subduction channel are neither prescribed nor assumed to represent a steady state. Instead, the system is free to evolve, starting from an imposed early stage of subduction, being controlled by the progressive modification of the thermal, petrological, and rheological structure of the subduction zone. In this evolution, upward migration of the aqueous fluid released from the subducting slab and progressive hydration of the mantle wedge play an important role. The results of our simulations are then compared to the record of some high-pressure metamorphic areas, including the Franciscan complex in California, the Yukon-Tanana terrane in northwestern Canada, and parts of the European Alps.

3. Methodological Background of the Simulation
3.1. Principal Equations [6] In this study, the numerical simulation is based on a set of differential equations describing the behavior of incompressible heat-conducting viscous media in the gravitational field. This set of equations comprises the Stokes equation of motion
@ sij =@ xj ? @ P=@ xi À rgi ; ð1Þ

the continuity equation
divðvÞ ? 0; ð2Þ ð3Þ
À1

and the heat transport equation
@ T =@ t þ vi @ T =@ xi ? ?@ ðk @ T =@ xi Þ=@ xi þ @ H =@ t =ðrCpÞ;

where x denotes the coordinates in m, v velocity in m Á s ; t _ time in s; sij =2heij the viscous deviatoric stress tensor in Pa; _ eij = 1/2(@ vi/@ xj + @ vj/@ xi) the strain rate tensors in sÀ1; P the pressure in Pa; T the temperature in K; h the effective viscosity in Pa Á s; r the density in kg Á mÀ3; g denotes the acceleration within the gravity field, 9.81 m Á sÀ2; k is the heat conductivity coefficient in W Á mÀ1 Á KÀ1; Cp is the isobaric heat capacity in J Á kgÀ1 Á KÀ1; and @ H/@ t is the radiogenic heat production in W Á mÀ3. Assuming a low shear stress in the subduction zones (see for instance discussion by Peacock [1990b, 1996], which is supported by the record of exhumed rocks combined with experimental results [Stockhert and Å Renner, 1998; Stockhert, 2002], the contribution of shear Å heating to the thermal budget of subduction zones is taken to be small and thus neglected. Also, we do not consider adiabatic effects and the latent heat of phase transformations. 3.2. Deformation Mechanisms and Rheology 3.2.1. Deformation in the Brittle Regime [7] The brittle deformation of a material, which is assumed to occur by frictional sliding on preexisting fracture planes of appropriate orientation, is described by the


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X-3

Mohr-Coulomb criterion [Brace and Kohlestedt, 1980; Ranalli, 1995]
t
max

? ðM1 sn þ M2 Þð1 À lÞ

ð4Þ

where tmax is the critical shear stress at which frictional sliding takes place given in Pa; sn is the normal stress on the fracture plane in Pa; l = Pfluid/Plith is the pore fluid pressure coefficient, i.e., the ratio between pore fluid pressure, Pfluid, and lithostatic pressure, Plith; M1 and M2 MPa are empirical constants (M1 = 0.85, M2 = 0 when sn < 200 MPa and M1 = 0.6, M2 = 60 MPa when sn > 200 MPa [Brace and Kohlestedt, 1980]. [8] For porous or fractured media containing a fluid phase, brittle strength is controlled by pore fluid pressure. For the upper crust, a hydrostatic pore fluid pressure gradient with a pore pressure coefficient l = 0.4 is generally accepted [e.g., Sibson, 1990]. Hydrocarbon exploration wells have shown that in sedimentary basins the transition from a hydrostatic to a near-lithostatic pore pressure gradient generally occurs at a depth of about 3 to 5 km [e.g., Sibson, 1990]. In contrast, the KTB drill hole has shown that in metamorphic basement rocks a hydrostatic pore pressure gradient can reach down to at least 9 km depth and a temperature of 265œC [Huenges et al., 1997; Grawinkel and Stockhert, 1997]. For simplicity, in the present Å simulation we assume a continuous transition from a hydrostatic pore fluid pressure (l = 0.4) at the surface to either a lithostatic pore fluid pressure (l = 0.99) in rocks for which a fluid-filled pore space is predicted, or to dry conditions (l = 0) in rocks devoid of a pore fluid, at a depth of 10 km. Intermediate pore fluid pressures are thus assumed at shallow depths <10 km, with an effective pore fluid pressure calculated as follows
l ? ?að10000 À zÞ þ bz=10000; ð5Þ

_ __ where eII ? 1=2eij eij is the second invariant of the strain rate tensor, with dimension sÀ1, F is a dimensionless coefficient dependent on the type of experiments on which the flow law is based (e.g., F = 2(1Àn)/n/3(1+n)/2n for triaxial compression and F = 2(1À2n)/n for simple shear). Dislocation creep is the deformation mechanism that poses an upper limit to strength for a given strain rate and temperature. Deformation in subduction zones, however, probably takes place by diffusion creep or dissolution precipitation creep in the presence of a fluid phase at a significantly lower stress for a given strain rate [e.g., Stockhert and Renner, 1998, Å Stockhert, 2002]. Å 3.2.3. Deformation by Diffusion Creep [10] Experimentally calibrated flow laws appropriate for diffusion creep and applicable to natural polyphase rocks are not available. However, Newtonian behavior and an inverse proportionality of strain rate and grain size d to the second or third power are predicted [e.g., Evans and Kohlstedt, 1995]. A generalized flow law then has the form _ g ? AN d
Àm

s expðÀE =RT Þ;

ð8Þ
À1

where z is the depth in m, and a = 0.4 is the pore fluid pressure factor in the hydrostatic regime, as a limiting case, while b is set either as 0.99 or as zero, assuming a quasi lithostatic pore fluid pressure [e.g., Etheridge et al., 1984] or a dry rock devoid of pore fluid at >10 km depth, respectively. For a wet rock, the pore pressure coefficient l increases from 0.4 at the surface to 0.99 at z = 10 km, implying continuous reduction of the brittle strength, while in the case of a dry rock the apparent strength increases with depth, with l approaching zero at z = 10 km. 3.2.2. Deformation by Dislocation Creep [9] For materials undergoing deformation by dislocation creep, the strength can be predicted using experimental flow laws of the power law type [e.g., Evans and Kohlstedt, 1995; Ranalli, 1995]
_ g ? AD sn expðÀE =RT Þ; ð6Þ _ where s is differential stress given in Pa, g is strain rate in sÀ1, E is the activation energy in kJ Á molÀ1, AD is a material constant with the dimension PaÀn Á sÀ1, n is the stress exponent, and R is the gas constant. A creep viscosity, hcreep, dependent on stress and temperature, is then defined in term of invariants by [Willett, 1999] h
creep

where AN is a material constant with the dimension Pa Á mm Á sÀ1; d is the grain size in m; and m is the grain size exponent (m = 2 or 3, depending on the precise mechanism and rate controlling step; see Evans and Kohlstedt [1995]). [11] For rocks that are expected to undergo deformation by diffusion creep, we assume a constant creep viscosity hcreep, independent of temperature. In the absence of appropriate calibrations of a flow law for polyphase material, this simplification appears to be justified by the fact that the viscosity in diffusion creep decreases with increasing temperature and increases with increasing grain size to the power of 2 or 3. In turn, grain size in metamorphic rocks tends to increase with increasing temperature, so that both effects counteract. For instance, assuming that the dependence of grain size on temperature (Figure 1a) can be approximated by the empirical relation
d ? expðÀ1:04 À 6001=T Þ; ð9Þ

predicting a characteristic grain size of about 0.01 mm at 300œC, of about 0.25 mm at 550œC, and of about 1 mm at 750œC, the viscosity would remain within one order of magnitude between 300 and 750œC for an activation energy of 150 kJ Á molÀ1 and a grain size exponent of 3 (Figure 1b). 3.2.4. Definition of an Effective Rheology [12] Following Schott and Schmeling [1998], the creep rheology is combined with a quasi-brittle rheology to yield an effective rheology. For this purpose the Mohr-Coulomb law (equation (4)) is simplified to the yield stress, syield criterion and implemented by a ``Mohr-Coulomb-viscosity,'' hMC as follows
h
MC

_ ? syield =ð4eII Þ
lith

1=2

;

syield ? ðM1 P

þ M2 Þð1 À lÞ:

ð10Þ

The total effective viscosity, h, is then defined by the following criterion
h?h
creep

_ when 2ðeII Þ _ when 2ðeII Þ

1=2

h

creep

< syield ; > syield

ð11aÞ ð11bÞ

_ ? ðeII Þ

ð1ÀnÞ=2n

F ðA D Þ

À1=n

expð E=nRT Þ;

ð7Þ

h?h

1=2

MC

h

creep


X-4

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 1. Empirical relation between effective grain size of metamorphic rocks and temperature (a), and effect on viscosity in grain-size sensitive diffusion creep (b), calculated for activation energies of 100, 150, and 200 kJ Á molÀ1. Note that the effect of temperature and temperature-dependent grain size cancel each other for an activation energy of 150 kJ Á molÀ1 over a wide temperature range. See text for discussion. 3.2.5. Mechanical Properties Assigned to the Different Materials [13] High-pressure metamorphic rocks exhumed from subduction zones reveal little evidence of deformation of their predominating minerals at depth by dislocation creep [Stockhert and Renner, 1998; Renner et al., 2001; Mauler et Å al., 2001; Sto khert, 2002]. Moreover, many high- and Åc ultrahigh-pressure metamorphic rocks have remained nearly undeformed. This requires that deformation during subduction and exhumation was localized along weak shear zones. Dissolution precipitation creep and fluid assisted granular flow are suspected to be the dominant deformation mechanisms in these shear zones, thus suggesting a Newtonian behavior. Combining an assumption on the cumulative width of the shear zones, typical plate velocities, and an upper bound to stress provided by the flow laws for dislocation creep of rock-forming minerals such as quartz and coesite [Renner et al., 2001; Stockhert, 2002], maxÅ imum bulk viscosities on the order of 1018 to 1020 Pa Á s are predicted. Apart from these rough estimates, little is known about the rheology of materials in subduction zones. Therefore, in our simulation, a Newtonian rheology is assigned to subducted sediment, the upper hydrated portion of the

basaltic oceanic crust, and to serpentinized mantle, while a power law rheology is used for the lower gabbroic portion of the oceanic crust, and for unhydrated mantle material. [14] As appropriate flow laws allowing more realistic predictions are not available, the Newtonian rheology is implemented by using a constant (Figure 1) effective viscosity hcreep with absolute values of 1018, 1019, or 1020 Pa Á s, used in different numerical experiments. This simplification is based on the assumption that deformation takes place by grain size and temperature dependent diffusion creep (equation (8)), with an activation energy of c. 150 kJ Á molÀ1 (which may be taken as a reasonable assumption) [e.g., Ranalli, 1995] and using the simple empirical dependence of grain size on temperature proposed above (equation (9)). The absolute value of the viscosity is thus considered as a variable parameter, which affects the structural evolution in the numerical simulation. [15] The power law flow law parameters are taken, for simplicity, from the compilation of Ranalli [1995], as the present simulations are considered as robust with respect to the choice of the power law parameters. [16] The rheology of the lower gabbroic layer of the oceanic crust is taken to be represented by a flow law for dislocation creep of plagioclase with composition An75, with E = 238 kJ Á molÀ1, n = 3.2, and logAD = À3.5 (AD given in MpaÀn Á sÀ1), as quoted by Ranalli [1995]. For this material a high pore fluid pressure is assumed (b = 0.99 in equation (5)). [17] Hydrated mantle rocks are expected to be much weaker than dry peridotite, caused by the presence of serpentine minerals and, at higher temperatures, of amphiboles. While olivine rheology is constrained by laboratory tests [e.g., Chopra and Paterson, 1981; Mei and Kohlstedt, 2000], with a pronounced effect of chemical environment [e.g., Bai et al., 1991; Mei and Kohlstedt, 2000], the mechanical properties of serpentine minerals and amphiboles are insufficiently constrained for high temperatures and low strain rates. Though serpentinized peridotite is generally assumed to be a very weak material [e.g., Hermann et al., 2000; Guillot et al., 2000], the available laboratory experiments in the brittle field [e.g., Escartin and Hirth, 1997; Escartin et al., 2001; Moore et al., 1997] are not necessarily appropriate to predict rheological properties of serpentinized mantle at high pressure and temperature, and at slow geological strain rates in a mantle wedge. More importantly, the mechanical properties of polyphase materials in the presence of a hydrous fluid cannot be predicted from presently available conventional laboratory tests. Nonetheless, the structural and microstructural record of natural rocks reveals a pronounced viscosity contrast between single phase and polyphase materials undergoing deformation by dissolution precipitation creep [Schwarz and Stockhert, 1996; Stockhert, 2002], with a peculiar role Å Å of sheet silicates enhancing dissolution. The underlying mechanisms are not yet fully explored, though the specific role of interphase boundaries as sites of enhanced dissolution rates compared to grain boundaries has been verified in experiments by, e.g., Hickman and Evans [1995] and Bos and Spiers [2000]. The importance of interfaces between


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X-5

unlike minerals in diffusion creep has also been pointed out in a theoretical model by Wheeler [1992]. Based on these considerations, significant weakening of the mantle can be expected due to replacement of olivine and pyroxene by serpentine minerals in the presence of an aqueous pore fluid. Here we use the antigorite stability field [Schmidt and Poli, 1998] to define the boundary between a weak serpentinized mantle and a strong serpentine-free mantle. [18] The rheology of unhydrated mantle is represented by a flow law for dislocation creep of dry olivine, with E = 532 kJ Á molÀ1, n = 3.5, and logAD = 4.4 (AD given in MpaÀn Á sÀ1) [Ranalli, 1995]. Brittle strength is high, assuming the absence of a free pore fluid (b = 0 in equation (5)). [19 ] T he rheology of th e partially hydrated mantle beyond the antigorite stability field is represented by a flow law for dislocation creep of wet olivine, with E = 470 kJ Á molÀ1, n = 4, and logAD = 3.3 (AD given in MpaÀn Á sÀ1) [Ranalli, 1995], assuming a high pore fluid pressure (b = 0.99 in equation (5)). [20] A Newtonian rheology with a constant viscosity of either 1018, 1019, or 1020 Pa Á s, independent of temperature (see above), is chosen in the different numerical experiments to represent the mechanical behavior of the serpentinized mantle within the antigorite stability field. 3.3. Numerical Techniques [21] The finite difference method on a regular rectangular 140 Â 100 non-deforming grid is used to solve equations (1) - (3). Pressures are calculated for the geometrical centers of cells, while temperatures and velocities are computed for the nodes (half-staggered grid). The marker technique [e.g., Weinberg and Schmeling, 1992] is used to move the field parameters (e.g., density or viscosity) with sharp discontinuities through a discrete mesh. Marker points containing material property information are initially distributed on the dense (700 Â 500) regular rectangular marker grid with small ( 1=2 of marker grid distance) random displacement. The markers are moved through the mesh according to the velocity field using the forth order Runge-Kutta method. [22] A newly developed code [Gerya et al., 2000] allows modeling of the P-T-time path for each given marker by using the calculated temperature and pressure fields at each t = const. This provides the possibility of investigating the shape of such a path and the field of P-T-time paths in the simulation. A detailed account of the method and tests of the accuracy of the numerical experiments is provided by Gerya et al. [2000].

[24] The thermal structure at the onset of subduction was modeled by Cloos [1982, 1985] and Peacock [1987a]. These studies demonstrate the evolution of the pattern of isotherms. At the onset of subduction, heat is rapidly withdrawn from the hanging wall of the subduction zone and transferred into the cool descending slab. With ongoing plate convergence, the hanging wall gradually becomes cooler until a steady state thermal configuration is reached. Accordingly, at an early stage, the subducted material follows a trajectory characterized by a higher T/P ratio, compared to the later stage steady state situation. In turn, the P-T path for hanging wall material in an immature system is characterized by isobaric cooling due to the heat withdrawn from the hanging wall. The exhumation path of such earlyaccreted material [e.g., Parkinson, 1996] is then characterized by lower T/P ratios compared to its burial path. The counterclockwise P-T trajectories found for the high-pressure metamorphic rocks of the Franciscan complex in California [Krogh et al., 1994] and for the Faro eclogite of the Yukon-Tanana terrane in the western Canadian Cordillera [Perchuk et al., 1999; Perchuk and Philippot, 2000b] are consistent with this scenario, and are further addressed in the discussion. 4.2. Initial and Boundary Conditions of the Simulations [25] The initial conditions chosen for our numerical experiments are displayed in Figure 2a. The subduction zone is prescribed by a weak layer [e.g., Kincaid and Sacks, 1997; Schmeling et al., 1999; van Hunen et al., 2000] of uniform thickness of 4 km and with an appropriate inclination. In view of the possible role of water in the initiation of subduction [Regenauer-Lieb et al., 2001] this weak layer is prescribed as a brittle fault within mantle rocks, characterized by the chosen wet olivine rheological parameters (see above), and a high pore fluid pressure is assumed, i.e., b = 0.99 in equation (5). During the progress of subduction, this pre-defined fault zone is substituted by weak subducted crustal rocks and serpentinized mantle, implying decoupling along the plate interface [e.g., Kincaid and Sacks, 1997; Schmeling et al., 1999; van Hunen et al., 2000]. This weak layer precludes significant deformation of the lithospheric portion of the mantle wedge, as observed in other numerical experiments [e.g., Pysklywec et al., 2000; Pysklywec, 2001]. [26] The initial temperature field (Figure 2a) is defined by an oceanic geotherm [Turcotte and Schubert, 1982] for a given lithospheric age, which varies between 20 and 80 Myr in the different numerical experiments. [27] Figure 2b shows the boundary conditions used in our model. These boundary conditions correspond to the corner flow model [e.g., Cloos, 1982; Guinchi and Ricard, 1999], modified to account for asthenospheric mantle flow at temperatures exceeding 1000œC in the mantle wedge [e.g., Davies and Stevenson, 1992; Kincaid and Sacks, 1997; van Hunen et al., 2000]. Following Davies and Stevenson [1992] we used the condition of zero deviatoric normal stress on two intersecting planes in order to obtain the velocity field in the hanging wall asthenosphere similar to the corner flow solution for an infinite wedge [e.g., Batchelor, 1967; Turcotte and Schubert, 1982]. The following

4. Model Design and Limitations
4.1. Constraints on the Early Stage of Subduction [23] The early stages in the evolution of a new subduction zone are poorly understood, as there appear to be no wellconstrained natural situations allowing a direct geophysical analysis. Insight into the early stages of subduction is thus restricted to either appropriate numerical experiments or the record of rocks exhumed from ancient subduction zones.


X-6

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

proceeding in the mantle wedge [e.g., Peacock, 1987b]. Therefore in our numerical experiments the hydration of the mantle wedge is modeled to proceed at a time- and spacedependent rate. The rate is considered to be a function of the distance from the trench, and of the total amount of fluid released from the plate [e.g., Peacock, 1987b, 1990b]. To account for the hydration process in a deforming medium we describe the vertical displacement of the hydration front with respect to the upper interface of the subducting slab (in Eulerian coordinates) by the following transport equation
@ zh =@ t ? vz À vx @ zh =@ x À vh ; ð12Þ

Figure 2. Definition of (a) initial and (b) boundary conditions used in the numerical experiments. conditions were defined for the lower right corner of the model: (1) zero deviatoric normal stress on the vertical plane at the right boundary and (2) zero deviatoric normal stress on the plane normal to the surface of the subducting plate (a plane inclined at a + 90œ; see Figure 2b) at the lower boundary. To account for the continuity of the temperature field across the lower boundary of the model we used an x- and time-dependent heat flow, q, corresponding to the condition @ 2q/@ z2 = 0. Test calculations have shown that this boundary condition yields a close approximation to the analytical solution for the relaxation of an infinite thermal profile in a static regime [e.g., Turcotte and Schubert, 1982]. In the dynamic regime of our numerical experiments the massive outward flow of matter means that the changes in the lower boundary condition for temperature (e.g., the use of the condition q = const [Guinchi and Ricard, 1999] do not significantly affect the temperature field inside the model space. 4.3. Dehydration and Hydration Reactions [28] In our numerical experiments we used a simple model (Figure 3) to describe the dehydration of the subducted slab and the related hydration of the mantle wedge, following the schemes proposed by Peacock [1993, 1996] and Schmidt and Poli [1998]. The availability of H2O is suggested to control the progress of the hydration reaction

where zh is the vertical position of the hydration front as a function of the horizontal distance, x, measured from the trench; vh is the hydration rate; vz and vx are the vertical and horizontal components of the material velocity vector at the hydration front. At each time step equation (12) was solved numerically to calculate the corresponding vertical displacement of the hydration front with respect to the upper surface of subducting plate. [29] In a situation where the availability of water controls the progress of the mantle hydration [e.g., Peacock, 1987b], spatial changes in the hydration rate should mainly depend on spatial changes in the rate of fluid release along the surface of the subducting plate. Following the arguments of Schmidt and Poli [1998], who favor a continuous rather than a single phase dehydration model for the subducting lithosphere, we assume that the rate of fluid release along the slab [e.g., Kerrick and Connolly, 2001] and, consequently, the hydration rate can be roughly approximated as a linear function of the horizontal distance, x, measured from the trench (Figure 3)
vh =vs ? A?1 À BtgðaÞx=zlim when tgðaÞx < zlim ;
lim

ð13aÞ ð13bÞ

vh =vs ? 0 when tgðaÞx > z

;

where a is the angle of inclination of the slab; vs is the subduction rate; zlim is the limiting depth below which fluid release from the subducting plate is taken to be negligible

Figure 3. Model for the hydration of the mantle wedge used in the numerical experiments. See text for details.


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS Table 1. Parameters of Numerical Models Studied
Model A B C D E F G H I J K L M Rate of Convergency, km/Myr 30 30 10 200 30 30 30 30 30 30 150 30 30 Subduction Angle, a 45 45 45 45 45 45 45 45 30 60 60 45 45
a

X-7

Initial Age of the Lithosphere, Myr 40 40 40 40 40 40 20 80 40 40 40 40 40

Viscosity Variations, Pa Á s 1019 -1024 1019 -1024 1019 -1024 1018 -1023 1019 -1024 1019 -1024 1019 -1024 1019À1024 1019À1024 1019À1024 1018À1023 1018À1024 1020À1024 A A A A A A A A A A A A A = = = = = = = = = = = = =

Hydration Progress: Parameters of Equations (13a), (14) 0.1, B = 1.0 (MH2O = 0.05, B = 0 (MH2O = 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.03 À 0.12, B = 1.0 0.09 À 0.35, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.07 À 0.25, B = 1.0 0.2 - 0.8 Â 106 kg/m2) 0.2 À 0.8 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.25 Â 106 kg/m2) (MH2O = 0.7 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2) (MH2O = 0.5 Â 106 kg/m2)

a Material properties in all studied models are taken as follows: heat conductivity coefficient,k (W mÀ1 KÀ1): 2 for all types of rocks; isobaric heat capacity, Cp (J kgÀ1 KÀ1): 1000 for all types of rocks; radiogenic heat generation, @ H/@ t (10À6 W mÀ3): sediments, 2; basaltic and gabbroic crust, 0.25; mantle, 0.022; density (constant, P-T independent), r (kg mÀ3): sediments, 2900; basaltic and gabbroic crust, 3100; nonserpentinized mantle, 3300; serpentinized mantle, 3000; rheological properties: sediments, basaltic crust and serpentinized mantle - newtonian creep (for all studied models an effective viscosity hcreep for the newtonian creep is taken to be equal to the lower limit of viscosity variations), gabbroic crust, An75 creep [Ranalli, 1995] with brittle strength calculated for high pore fluid pressure (b = 0.99 in equation (5)); unhydrated mantle - dry olivine creep [Ranalli, 1995] with brittle strength calculated for dry conditions (b = 0 in equation (5)), strong (nonserpentinized) hydrated mantle - wet olivine creep [Ranalli, 1995] with brittle strength calculated for high pore fluid pressure (b = 0.99 in equation (5)).

[e.g., Peacock, 1990b; Phillipot, 1993; Schmidt and Poli, 1998; Scambelluri and Philippot, 2001]. The limiting depth depends on the amount of fluid stored in the subducting lithosphere and on the progress of dehydration and partial melting, respectively progressive dissolution of material in supercritical fluids [e.g., Shen and Keppler, 1997; Bureau and Keppler, 1999]. The limiting depth has been proposed to vary between ca. 50 and >300 km as a function of the geotherm [e.g., Phillipot, 1993; Schmidt and Poli, 1998]. For simplicity, here the limiting depth is defined as the depth corresponding to the P-T conditions of the wet basalt solidus [Schmidt and Poli, 1998] within the subducted oceanic crust (Figure 3). It is assumed that beyond this depth the fluid phase in subducted oceanic crust is represented by a high density and viscosity partial melt [e.g., Peacock, 1990b] that does not produce hydration of the mantle wedge. In all models presented here, this depth is either greater or very close to the maximum depth of antigorite stability in the mantle wedge. This condition is valid for relatively warm subduction zones, where the P-T conditions of the wet basalt solidus may be reached in the oceanic crust at about 100 km depth [e.g., Schmidt and Poli, 1998]. [30] Parameter A in equation (10) may be interpreted as a non-dimensional intensity of the process of hydration of the mantle wedge. The non-dimensional parameter B may vary from 0 to 1, characterizing the decrease in hydration rate with depth. By assuming that the quantity of water released from the subducted lithosphere equals that consumed by hydration of the mantle wedge [e.g., Schmidt and Poli, 1998], parameter A in equation (13) becomes
A ? 2tgðaÞMH2O =?ð2 À BÞrw WH2O zlim ; ð14Þ

hydrated mantle, given as kg mÀ3; WH2O is the weight fraction of water in the hydrated mantle. Using the relations proposed by Schmidt and Poli [1998], we relate the weight fraction of water to the density of the hydrated mantle rw by the empirical equation
WH2O ? 0:07ð3300 À rw Þ=420: ð15Þ

[31] Taking into account the possible variations in the parameters of equation (14) (MH2O ranging between 0.25 Â 106 and 0.70 Â 106 kg/m2; zlim ranging between 70 Â 103 and 300 Â 103 m; rw ranging between 2880 and 3150 kg/m3 [Schmidt and Poli, 1998], the typical values of the parameter A range between 0.01 and 0.30. In our numerical experiments, the influence of these variations in A (as a function of MH2O for rw = 3000 kg/m3) and B on the rate and extent of the hydration of the mantle wedge is investigated as a function of position. Then, the effect on the particle paths and the thermal structure of the subduction zone is explored, and thus the effect on the P-T-time paths recorded by hypothetical high-pressure metamorphic rocks. [32] The densities of all rock types are taken as constant and P-T independent (Table 1, comments). The only change in density considered in our model is related to the transition between unhydrated (3300 kg/m3) and partly serpentinized (3000 kg/m3) mantle material.

5. Results of the Simulations
5.1. Evolution of the Subduction Zone [33] The structure and the circulation pattern of the developing subduction zone in our numerical simulations are depicted in Figures 4 and 5. The evolution can be subdivided into two different stages: (1) the early stage, without return flow from depths exceeding about 20 km (corresponding to the time slices between 1.8 Myr and 10.7

where MH2O denotes the total amount of water released from the subducting plate by the limiting depth, given in kg mÀ2; a is the subduction angle; rw is the density of the


X-8

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 4. Development of a wedge-shaped circulation pattern of material during ongoing subduction. The initial and boundary conditions correspond to 140 Â 100 km model A described in Table 1. Left: Evolution of the temperature field (isotherms labeled in œC) and distribution of rock types (color code): 1 = sedimentary rocks; 2 = basaltic upper oceanic crust; 3 = gabbroic lower oceanic crust; 4 = unhydrated mantle; 5 = serpentinized mantle; 6 = hydrated mantle beyond stability field of antigorite. Right: Evolution of viscosity (color code) and velocity field (arrows). See color version of this figure at back of this issue. Myr), and (2) the mature stage, with intense return flow from greater depth, involving the hydrated portion of the mantle wedge (corresponding to the time slices between 10.7 and 32.8 Myr). [34] The transition between these two stages is controlled by the progress of mantle wedge hydration. In the early stage (time slices 1.8 to 10.7 Myr), a shallow (<20 km) circulation pattern develops within the accretionary wedge close to the trench. In the mature stage (time slices 10.7 to 24.6 Myr) an independent deep-seated circulation pattern has developed at 30 to 70 km depth in the hydrated portion of the mantle wedge, providing an effective mechanism for the exhumation of high-pressure metamorphic rocks. Still later (between 24.4 and 32.8 Myr) in the mature stage, further serpentinization of the mantle lithosphere of the upper plate results in coupling of both systems and in formation of a tectonic window with a high-pressure metamorphic melange. [35] A characteristic feature observed in all simulations is the widening of the subduction channel caused by progres-


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X-9

Figure 5. Development of a channel pattern of circulation of material during ongoing subduction. The initial and boundary conditions correspond to 140 Â 100 km Model B in Table 1. Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue. sive hydration of the mantle wedge, giving rise to a complex internal geometry, with a tectonic melange composed of hydrated fragments of the mantle wedge and subducted oceanic crust. The proportion of oceanic crust involved in circulation in the subduction channel strongly depends on the choice of the hydration rate as a function of position (Figure 3), which affects the shape and extent of the hydrated part of the mantle wedge. The highest proportion of incorporated oceanic crust is observed for a wedgeshaped channel (Figure 4), as a result of a marked decrease in the hydration rate with increasing horizontal distance from the trench (B = 1 in equation (13)). In contrast, simulations with a plane parallel shaped hydrated zone (Figure 5) as a result of a constant hydration rate (B = 0 in equation (13)), reveal comparatively little incorporation of oceanic crust. [36] In the majority of the simulations, the closure of the subduction channel at depth is located at the level where forced asthenospheric flow in the mantle wedge becomes effective (Figures 4 and 5). This zone is characterized by steep horizontal gradients in temperature, viscosity, and stress state. In the case of the planar channel the decreasing width of the hydrated zone between the top of the slab and the hydration front is controlled by the particle velocity field resulting from the interaction the two realms of viscous flow (Figure 5). The wedge-shaped geometry of the closure of


X - 10

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 6. Variations in geometry and thermal structure of the subduction zone predicted by the numerical models after 500 km of convergence, depending on the rate of subduction (a and b) and the characteristics of the hydration process (c and d), see text for details. The initial and boundary conditions correspond to models C (a), D (b), E (c) and F (d) described in Table 1. The inset separate sketches represent enlarged 110 Â 100 km areas of the original 140 Â 100 km models. Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue. the subduction channel causes return flow in the hydrated portion of mantle wedge. [37] The large-scale structure and the temperature field evolving in the numerical simulations are controlled by three types of material flow: subduction flow, asthenospheric flow beneath the overriding plate, and return flow in the hydrated portion of the mantle wedge. Interaction between these realms causes a complex evolution of the flow patterns and the thermal structure. The asthenospheric flow beneath the overriding plate and the return flow in the hydrated portion of the mantle wedge significantly affect the evolution of the temperature field (Figures 4 and 5). Compared to the models of Peacock [1990a, 1996], our simulations reveal a more complex thermal structure and slower cooling of the mantle wedge, despite the inclination of the subduction zone and rate of convergence chosen being identical in both models. 5.2. Sensitivity to Important Parameters [38] Table 1 summarizes the parameters used in our set of numerical simulations. The results of the experiments after a convergence of 500 and 1000 km are shown in Figures 6, 7, and 8. Figures 6a and 6b show that, for the same total convergence, a convergence rate of 1 cm/year results in a deeper reaching hydration zone and a larger proportion of oceanic crust being incorporated into the subduction channel, compared to subduction with a convergence rate of 20 cm/year. On the other hand, variations in the hydration rate do not significantly affect the absolute amount of incorporated oceanic crust, although the ratio between oceanic crust and serpentinized mantle rocks in the tectonic melange is modified (compare Figures 6c and 6d). The depth of the lower limit of the hydrated zone depends on the temperature field and increases with time and progressive cooling of the mantle wedge (Figures 4 and 5). Accordingly, at the same total convergence, subduction of 20 Myr old lithosphere results in a shallower lower limit of the hydrated zone compared to subduction of 80 Myr old lithosphere (compare Figures 7a and 7b). The angle of inclination of the subducting plate affects the thermal structure of the mantle wedge and the position of the tectonic window, composed of exhumed high-pressure metamorphic tectonic melanges, relative to the trench (compare Figures 7c and 7d). Figure 8 illustrates the influence of the rheology of the serpentinized mantle on the internal structure of the tectonic melange. The lower the effective viscosity of the serpentinized mantle, the higher is the degree of disintegration and mixing of different rock types. The position of the tectonic window is also modified dependent on the viscosity. 5.3. P-T-Time Trajectories of Rocks [39] The P-T-time paths of high-pressure metamorphic rocks that returned into the upper crust in our numerical simulations are highly variable (Figure 9). This is related to both the shape of the trajectories of circulation of these rocks in the hydrated portion of the mantle wedge (Figures 9b and 9c), and to the evolution of the temperature field with progressive subduction (Figure 4). Remarkably, both clockwise (i.e., exhumation at higher temperatures compared to burial) and counterclockwise (i.e., exhumation at


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X - 11

Figure 7. Variations in the geometry and the thermal structure of the subduction zone after 500 km of convergence, depending on the initial age of the lithosphere (a and b) and on the inclination of the slab a (c and d). The initial and the boundary conditions correspond to models G (a), H (b), I (c) and J (d) specified in Table 1. The inset separate sketches represent enlarged 110/tg(a) Â 100 km areas of the original 140/tg(a) Â 100 km models. Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue. lower temperatures compared to burial) P-T paths can be found for exhumed high-pressure metamorphic fragments of oceanic crust (Figures 9 and 10). Counterclockwise P-T paths are characteristic of discontinuous circulation of oceanic material subducted in a very initial stage and accreted to the hanging wall at great depth (compare the open square markers in Figures 10 and 11). Rapid nearisobaric cooling of these rocks, with a rate of >30œC/Myr (Figures 11a and 11c), is then caused by the displacement of the isotherms to greater depth in the initial stage of subduction (Figure 10, time slices 6.4 to 15.3 Myr). At a later stage, these rocks can be set free by hydration of the mantle wedge (Figure 10, time slices 19.4 to 25.3 Myr) and return toward the surface. On the other hand, clockwise P-T paths are characteristic of the continuous circulation of fragments of oceanic crust within the hydrated mantle wedge (compare solid square markers in Figures 10 and 11). Exhumation of high-pressure metamorphic rocks revealing both types of P-T paths may be contemporaneous (Figure 10, time slices 19.4 to 25.3 Myr), and finally they may be juxtaposed at a shallow crustal level within the tectonic melange (Figures 9a and 10). Obviously, P-T trajectories for the exhumation of high-pressure rocks fall into a P-T field of stability of antigorite in the mantle wedge (Figure 9c). [40] Comparison of the burial and exhumation histories (Figure 11) of the high-pressure metamorphic rocks, for a constant moderate convergence rate of 30 km/Myr, reveals that 1. the maximum burial rates are 15 to 20 km/Myr, and are thus 2 to 3 times higher than the maximum exhumation rates of 5 to 8 km/Myr. This is consistent with analytical results obtained for the case where subduction flux in the channel equals return flux [Gerya and Stockhert, 2002]; Å 2. the maximum heating rates are 150 to 250œC/Myr, and are thus 5 to 8 times higher than the maximum cooling rates of 30 to 50œC/Myr. Both the burial and exhumation rates and the heating and cooling rates increase proportionally with increasing rate of convergence (Figures 12 and 13). Consequently, the quoted ratios between the maximum rates of burial and exhumation, respectively heating and cooling, do not vary significantly when the rate of convergence is changed (compare Figures 10 and 11, and Figures 12 and 13, respectively). Also, rapid return flow in the subduction channel is favored by a Newtonian rheology of the subduction channel media [Gerya and Stockhert, 2002]. When the subduction flux in Å the channel either equals or exceeds the return flux, as realized in our simulations, the maximum exhumation rates predicted for a Newtonian rheology are at least twice as high as those predicted for a power law rheology.

6. Comparison of the Predicted Length and Timescales With Two Natural Examples
6.1. Franciscan Complex [41] The Franciscan complex in California developed along the active continental margin between the North American plate and some subducted precursor of the Pacific plate in Mesozoic times. It is composed of several belts of broken formation or melange, with a matrix of fine-grained


X - 12

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

structural record of the blocks, which are weakly deformed [Cloos, 1986]. Our simulation implies that the incorporation of the higher-grade blocks into the low-grade metasedimentary matrix could be the result of interaction between the shallow circulating system representing the accretionary wedge, and the deep circulation in the subduction channel. The age of metamorphism of the higher-grade blocks has been generally found to be in the range of 165 to 140 Ma [e.g., Coleman and Lanphere, 1971; Suppe and Armstrong, 1972; Mattinson, 1986], while most apparent ages of lowgrade metamorphism of the metasedimentary matrix are between 130 and 120 Ma [e.g., McDowell et al., 1984]. This relative timing reconciles with the results of our simulations, with interaction between the shallow system and the material exhumed from the deep-seated subduction channel after ca. 30 Myr of evolution of the subduction system (Figures 4 and 10). 6.2. Zermatt Saas Zone, Western Alps [42] The Zermatt-Saas Zone in the western Alps is composed of slices of oceanic crust and mantle derived from the Jurassic Piemont-Ligurian ocean subducted prior to

Figure 8. Variations in the geometry of the subduction zone after 1000 km of convergence, depending on the viscosity of serpentinized mantle. The initial and boundary conditions correspond to models L (a), and M (b) specified in Table 1. The separate sketches represent enlarged 110 Â 100 km areas of the original 140 Â 100 km models. Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue.

clastic sedimentary rocks and greywackes, that in places have recorded high-pressure metamorphism at low temperatures [e.g., Ernst, 1993; Terabayashi et al., 1996], with minor cherts and carbonate rocks. Blocks of serpentinites and higher-grade high-pressure metamorphic rocks, mainly amphibolites, eclogites, and blueschists, are locally embedded in the matrix. The size of these blocks and slices is generally small [Coleman and Lanphere, 1971], between a few and some hundred meters, and their grade of metamorphism is highly variable [e.g., Moore, 1984; Moore and Blake, 1989; Krogh et al., 1994; Oh and Liou, 1990]. Notably, some of the amphibolites and eclogites reveal counterclockwise P-T paths [e.g., Krogh et al., 1994; Wakabayashi, 1990]. The small size of the higher-grade metamorphic blocks, invariably derived from oceanic lithosphere, and the variability of the recorded P-T paths has been taken to suggest an origin in a subduction channel with forced flow [Cloos, 1982]. Notably, most of the highergrade metamorphic blocks are reported to show evidence of incorporation into a serpentinitic matrix at an earlier stage of their history [Cloos, 1986]. The low viscosity of the material in the subduction channel is consistent with the

Figure 9. Particle trajectories and P-T paths calculated for the model with a wedge-shaped circulation pattern. The P-T paths represent the history of exhumed high-pressure metamorphic rocks. The initial and boundary conditions correspond to Model A specified in Table 1. The sketch in (a) represents the original 140 Â 100 km model geometry for the instant of ``sampling'' of the particle trajectories (b) and the corresponding P-T paths (c). Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue.


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X - 13

Figure 10. Contrasting P-T evolution of two fragments of oceanic crust predicted for subduction with a moderate rate of 30 km/Myr. The initial and boundary conditions correspond to model A specified in Table 1. Separate sketches represent enlarged 74 Â 71 km areas of the original 140 Â 100 km model. Colors and symbols correspond to those given in Figure 4. See color version of this figure at back of this issue. collision between the Adriatic and the European continental margins [Dal Piaz and Ernst, 1978; Dal Piaz et al., 1993; Polino et al., 1990]. It is exposed over about 20 km normal to strike, forming a composite nappe with a thickness of between less than 1 and about 2 km. It comprises slices of eclogitized basalt and gabbro, and extensive serpentinized peridotite. Within the Zermatt-Saas Zone, the recorded conditions of eclogite facies metamorphism are variable,


X - 14

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

conditions [Cartwright and Barnicoat, 2002]. Notably, the coesite-bearing eclogites exhumed from >80 km show no record of deformation between the stage of maximum burial and exhumation to less than 40 km [van der Klauw et al., 1997], which again reconciles with transport in a lowviscosity subduction channel.

7. Discussion and Conclusions
[43] Our numerical experiments indicate that one feasible model for the exhumation of high- and ultrahigh-pressure metamorphic rocks involves forced flow in an either wedgeshaped or planar subduction channel, with geometry progressively evolving due to hydration of the mantle wedge (Figures 4 and 5). The role of serpentinized mantle has previously been emphasized, based on field relations, by, e.g., Blake et al. [1995], Hermann et al. [2000], Guillot et al. [2000, 2001], Kasahara et al. [2001] and Bostock et al. [2002]. Also, Schwartz et al. [2001] recently discussed the possible role of serpentinites as a transport medium for the exhumation of eclogitic rocks, based on calculations assuming buoyancy-driven flow of serpentinized peridotites relative to a dense anhydrous mantle wedge. In our model, however, the buoyancy forces related to the density contrast between non- serpentinized and serpentinized mantle (Dr = 300 kg/m3) do not notably affect the flow pattern, as the unhydrated peridotites on top of the serpentinized portion of the channel reveal a very high effective viscosity ()1021 Pa s, Figures 4 and 5). Therefore, buoyancy-driven escape of serpentinites from the subduction channel does not play a significant role compared to displacement of material within the channel driven by forced flow. [44] The simulations imply that the volume ratio between high-pressure metamorphic crustal rocks and serpentinized mantle material depends on the progress of hydration in the mantle wedge which controls the incorporation of subducted crust into the subduction channel (e.g., Figures 6c and 6d). Unfortunately, the restrictions inherent in simple models for the complex process of the hydration of the mantle wedge [e.g., Peacock, 1987b; Phillipot, 1993; Bebout and Barton, 1993; Manning, 1996; Guillot et al., 2001; Scambelluri and Philippot, 2001] limit the applicability of our numerical experiments for understanding the details of the hydration process. In particular, the insufficient understanding of the transport and mechanical properties of serpentinites at high temperatures and great depth, and the lack of appropriate kinetic data required to predict the rates of hydration and dehydration reactions, hamper the development of more sophisticated physical models concerning the release, transport, and consumption of aqueous fluids in subduction zones. [45] Furthermore, our present numerical experiments do not account for erosion that may significantly affect the displacement and exhumation of high-pressure metamorphic rocks, especially at the relatively shallow (<30 km) level of the subduction zone [e.g., Guinchi and Ricard, 1999; Beaumont et al., 1999; Ellis et al., 1999]. [46] The rates of change of P and T predicted in our simulations can be compared with the record of natural

Figure 11. Temperature (a), pressure (b), heating respectively cooling rate (c), and burial respectively exhumation rate (d) calculated for two markers displayed in Figure 10.

scattering between ca. 450œC at 1.2 GPa, and 600œC at 2.8 GPa [Martin and Tartarotti, 1986; Barnicoat and Fry, 1986; Reinecke, 1991; van der Klauw et al., 1997; Cartwright and Barnicoat, 2002]. Burial to >80 km depth is demonstrated by the occurrence of coesite in both oceanic metasediments and eclogites at Lago Cignana [Reinecke, 1991, 1998; van der Klauw et al., 1997], but has so far not been identified elsewhere. Exhumation of this UHP metamorphic oceanic slice took place within a few million years [Amato et al., 1999]. The considerable spread in the recorded P-T conditions, and their rather unsystematic regional distribution within the unit, indicate that the Zermatt-Saas Zone does not represent a coherent oceanic slab, but is composed of slices w ith contrasting burial and exhumat ion paths, finall y agglomerated at a depth of less than about 40 km, reconciling with the results of our simulations (e.g., Figures 8b and 9). Again, identified shear zones within the Zermatt-Saas Zone appear to be invariably formed at greenschist facies


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X - 15

Figure 12. P-T evolution of two fragments of oceanic crust during subduction and exhumation. The initial and boundary conditions correspond to Model K specified in Table 1. The separate sketches represent enlarged 43 Â 71 areas of the original 81 Â 100 km model. Colors and symbols correspond to those given in Figures 4 and 9, respectively. See color version of this figure at back of this issue.

rocks. For example, the eclogites exposed near Faro in the Yukon-Tanana terrane of the western Canadian Cordillera [Erdmer et al., 1998; Philippot et al., 2001] experienced a prograde evolution from about 520œC and !1.1 GPa to

690œC and !1.5 GPa, followed by near isobaric cooling down to about 540œC [Perchuk et al., 1999]. The counterclockwise P-T path was interpreted to reflect the tectonic evolution during the earliest stage of subduction. Based on


X - 16

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 13. Temperature (a), pressure (b), heating respectively cooling rate (c), and burial respectively exhumation rate (d) calculated for two markers displayed in Figure 12. diffusion modeling on zoned garnet, Perchuk et al. [1999] have derived rates of burial and heating of 70 km/Myr and 950œC/Myr, respectively, and a minimum cooling rate of 750œC/Myr. These very high rates of P and T change deduced from the diffusion models are supported by the results of Lu-Hf and Ar-Ar dating on garnet and phengite. They compare well to the results of the numerical experiment for the case of rapid subduction (Figures 12 and 13) with a rate of 150 km/Myr and a relatively steep inclination of 60œ (Model K in Table 1), where burial rates of 50 to 110 km/Myr at T > 500œC and heating rates of 500 to 1300œC/Myr at T > 500œC are attained. Both the shape of the simulated P-T path and the absolute values of pressure and temperature for the different stages are consistent with the natural record. On the other hand, the cooling rates of 100 to 200œC/Myr found in the simulations for the stage of near isobaric cooling (Figure 13c) are lower, but still comparable within the limits of uncertainty to the rates

derived from the preservation of compositional zoning of garnet. [47] The flow patterns and structures predicted in our numerical simulations are controlled by the rheological properties assigned to the different materials. It should be noted that the assumption of a Newtonian rheology with a low effective viscosity is not restricted to the serpentinized mantle wedge, but in our simulations also applies to subducted sediment and hydrated upper oceanic crust. As Newtonian behavior and low effective viscosities due to deformation by some kind of dissolution precipitation creep or fluid-assisted granular flow are suspected to predominate in crustal materials buried in subduction zones in general [e.g., Stockhert, 2002], the predicted flow patterns may be Å of more general significance, possibly also applying to situations where continental crust is incorporated into the subduction channel. Striking similarities are evident when comparing the structures developed in the upper level of the subduction channel in our simulations (Figure 8b; also compare with nappe-like structures produced as the result of a forced flow within an orogenic wedge in numerical experiments by Allemand and Lardeaux [1997], and Guinchi and Ricard [1999]) with cross sections through the largely high-pressure metamorphic nappes of the western Alps [e.g., Debelmas et al., 1983; Polino et al., 1990; Dal Piaz et al., 1993]. It may be suspected that the internal structure of orogenic belts with exposed high-pressure and ultrahigh-pressure metamorphic units, commonly referred to as nappes, with a thickness on the order of 102 to 103 m and a horizontal width of exposure on the order of some 103 to 104 m, which feature different maximum P and T conditions, and contrasting P-T histories, may in fact partly represent material that has been extruded from a widening subduction channel to form an orogenic root beneath the forearc during ongoing subduction. Incorporated continental material can be derived from the active continental margin by subduction erosion [von Huene and Scholl, 1991] and need not always indicate collision or accretion of terranes or microcontinents from the subducted plate. Clearly, the final exhumation of such assemblies and exposure at the surface can be a consequence of later collision, with considerable post-high-pressure metamorphic deformation and obliteration of the earlier record. [48] Notwithstanding the difficulties of any geological application of results obtained by numerical simulations on the given length and time scales, the following conclusions appear to be robust and are easily reconciled with the natural record, though clearly are not necessarily unequivocal. 1. Burial and exhumation of high-pressure metamorphic rocks in subduction zones are likely affected by progressive hydration of the mantle wedge. This process controls the shape and internal circulation pattern of a subduction channel. Widening of the subduction channel due to hydration of the hanging wall mantle results in the onset of forced return flow in the subduction channel. This explains why the association of high- and/or ultrahighpressure metamorphic rocks with more or less hydrated mantle material is often characteristic for high-pressure metamorphic terranes.


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

X - 17

2. The shape of the P-T path, and the maximum P-T conditions achieved by an individual high-pressure metamorphic rock, depend on the specific trajectory of circulation in the subduction channel. Both clockwise and counterclockwise (i.e., with lower temperatures during exhumation compared to burial) P-T paths are possible for slices of oceanic crust that became involved in the circulation. Counterclockwise P-T paths are found for slices that are accreted to the hanging wall at an early stage of subduction, and set free by the progress of hydration and softening in a more evolved stage, returning towards the surface in a cooler environment. On the other hand, slices that were involved in continuous circulation, or that entered the subduction zone when a steady state thermal structure was already achieved, reveal exclusively clockwise trajectories. Our model also indicates that P-T trajectories for the exhumation of high-pressure rocks in subduction channel fall into a P-T field of stability of antigorite in the mantle wedge. 3. An array of diverse, though interrelated, P-T paths rather then a single P-T trajectory is expected to be

characteristic for subduction related metamorphic complexes. In fact, this appears to be what is commonly observed in orogenic belts. The characteristic size and shape of the units with an individual history depend on the rheology of the material in the subduction channel. Our simulations indicate that lower viscosities result in enhanced disintegration and mingling, and hence smaller characteristic length scales for coherent units and a marked contrasts between adjacent slices, a structure commonly termed melange, while higher viscosities favor the formation of extensive coherent slices, possibly representing the type of structures commonly referred to as nappes in high-pressure metamorphic terranes.

[49] Acknowledgments. This work was supported by RFBR grants 0005-64939 and 01-05-64585 and by Alexander von Humboldt Foundation Research Fellowships to T. V. Gerya and A. L. Perchuk, and by the Sonderforschungsbereich 526 at Ruhr-University, funded by the Deutsche Forschungsgemeinschaft. Jorg Renner is thanked for discussion and comÅ ments, and Stuart Thomson for checking the English. Constructive reviews by Chris Beaumont and Laurent Jolivet are greatly appreciated.

References
Allemand, P., and J.-M. Lardeaux, Strain partitioning and m et am o r phi sm in a d eformabl e o rogen i c wedge: Application to the Alpine belt, Tectonophysics, 280, 157 - 169, 1997. Amato, J. M., C. M. Johnson, L. P. Baumgartner, and B. L. Beard, Rapid exhumation of the Zermatt-Saas ophiolite deduced from high-precision Sm-Nd and Rb-Sr geochronology, Earth Planet. Sci. Lett., 171, 425 - 438, 1999. Bai, Q., S. J. Mackwell, and D. L. Kohlstedt, Hightemperature creep of olivine single crystals, 1, Mechanical results of buffered samples, J. Geophys. Res., 96, 2441 - 2464, 1991. Barnicoat, A. C., and N. Fry, High-pressure metamorphism of the Zermatt-Saas ophiolite zone, Switzerland, J. Geol. Soc. London, 143, 603 - 618, 1986. Batchelor, G. K., An Introduction to Fluid Dynamics, 615 pp., Cambridge Univ. Press, New York, 1967. Beaumont, C., S. Ellis, and O. A. Pfiffner, Dynamics of sediment subduction-accretion at convergent margins: Short-term modes, long-term deformation, and tectonic implications, J. Geophys. Res., 104, 17,573 - 17,601, 1999. Bebout, G. E., and M. D. Barton, Metosomatism during subduction: Products and possible paths in the Catalina Schists, California, Chem. Geol., 108, 61 - 92, 1993. Blake, M. C. Jr., D. E. Moore, and A. S. Jayko, The role of serpentinite melanges in the unroofing of UHPM rocks: An example from the Western Alps of Italy, in Ultrahigh Pressure Metamorphism, edited by R. G. Coleman and X. Wang, pp. 182 - 205, Cambridge Univ. Press, New York, 1995. Bos, B., and C. J. Spiers, Effect of phyllosilicates on fluid-assisted healing of gouge-bearing faults, Earth Planet. Sci. Lett., 184, 199 - 210, 2000. Bostock, M. G., R. D. Hyndman, S. Rondenay, and S. M. Peacock, An inverted continental Moho and serpentinization of the forearc mantle, Nature, 417, 536 - 538, 2002. Brace, W. F., and D. L. Kohlestedt, Limits on lithospheric stress imposed by laboratory experiments, J. Geophys. Res., 85, 6248 - 6252, 1980. Bureau, H., and H. Keppler, Complete miscibility between silicate melts and hydrose fluids in the upper mantle: Experimental evidence and geochemical implications, Earth Planet. Sci. Lett., 165, 187 - 196, 1999. Carswell, D. A., and R. Y. Zhang, Petrographic characteristics and metamorphic evolution of ultrahighpressure eclogites in plate-collision belts, Int. Geol. Rev., 41, 781 - 798, 1999. Cartwright, I., and A. C. Barnicoat, Petrology, geochronology, and tectonics of shear zones in the ZermattSaas and Combin zones of the Western Alps, J. Metamorph. Geol., 20, 263 - 281, 2002. Chopin, C., Coesite and pure pyrope in high-grade blueschists of the Western Alps: First record and some consequences, Contrib. Mineral. Petrol., 86, 107 - 118, 1984. Chopra, P. N., and M. S. Paterson, The experimental deformation of dunite, Tectonophysics, 78, 453 - 473, 1981. Cloos, M., Flow melanges: Numerical modelling and geologic constraints on their origin in the Fransiscan subduction complex, California, Geol. Soc. Am. Bull., 93, 330 - 345, 1982. Cloos, M., Thermal evolution of the convergent plate margins: Thermal modeling and reevaluation of isotopic Ar-ages for blueschists in the Franciscan complex of California, Tectonics, 4, 421 - 433, 1985. Cloos, M., Blueschists in the Franciscan complex of California: Petrotectonic constraints on uplift mechanisms, in Blueschists and Eclogites, edited by B. W. Evans and E. H. Brown, Mem. Geol. Soc. Am., 164, 77 - 93, 1986. Cloos, M., and R. L. Shreve, Subduction-channel model of prism accretion, melange formation, sediment subduction, and subduction erosion at convergent plate margins, 1, Background and description, Pure Appl. Geophys., 128, 455 - 500, 1988a. Cloos, M., and R. L. Shreve, Subduction-channel model of prism accretion, melange formation, sediment subduction, and subduction erosion at convergent plate margins, 2, Implications and discussion, Pure Appl. Geophys., 128, 501 - 545, 1988b. Coleman, R. G., and M. A. Lanphere, Distribution and age of high-grade blueschists, associated eclogites, and a mphi boli tes from O re gon an d C ali f or ni a, Geol. Soc. Am. Bull., 82, 2397 - 2412, 1971. Dal Piaz, G. V., and W. G. Ernst, Areal geology and petrology of eclogites and associated metabasites of the Piemonte ophiolite nappe, Breuil-St.Jacques ar ea, Italian Western Alps, Tectonophysics, 51, 99 - 126, 1978. Dal Piaz, G. V., G. Gosso, G. Pennacchioni, and M. L. Spalla, Geology of eclogites and related rocks in the Alps, in Italian Eclogites and Related Rocks, edited by L. Morten, pp. 17 - 58, Scritti e Doc. Accad. Naz. delle Sci. Detta dei XL, Roma, 1993. Davies, H. J., and D. J. Stevenson, Physical model of source region of subduction zone volcanics, J. Geophys. Res., 97, 2037 - 2070, 1992. Debelmas, J., A. Escher, a nd R. Tr u p y, Pr o fi l e s Åm through the western Alps, in Profiles of Orogenic Belts, Geodyn. Ser., vol. 10, edited by N. Rast and F. M. Delany, pp. 83 - 96, AGU, Washington, D.C., 1983. Ellis, S., C. Beaumont, and O. A. Pfiffner, Geodynamic models of crustal-scale episodic tectonic accretion and underplating in subduction zones, J. Geophys. Res., 104, 15,169 - 15,190, 1999. Erdmer, P., E. D. Ghent, D. A. Archibald, and M. Z. Stout, Paleozoic and Mesozoic high-pressure metamorphi sm at the margin of the ancestral North America in central Yukon, Geol. Soc. Am. Bull., 110, 615 - 629, 1998. Ernst, W. G., Mineral parageneses and plate tectonic settings of relatively high-pressure metamorphic belts, Fortschr. Mineral., 54, 192 - 222, 1977. Ernst, W. G., Metamorphism of Franciscan tectonostratigraphic assemblage, Pacheco Pass area, east-central Diablo Range, California Coast Ranges, Geol. Soc. Am. Bull., 105, 618 - 636, 1993. Ernst, W. G., Metamorphism, partial preservation, and exhumation of ultrahigh-pressure belts, Island Arc, 8, 125 - 153, 1999. Ern st, W. G., S ubdu ction , ultra high- pre ssur e me tam orp hism, a nd r egurgitation of buoya nt c rusta l slices--Implications for arcs and continental growth, Phys. Earth Planet. Inter., 127, 253 - 275, 2001. Ernst, W. G., and J. G. Liou, Overview of UHP Metamorphism and tectonics in well-studied collisional orogens, Int. Geol. Rev., 41, 477 - 493, 1999. Ernst, W. G., J. G. Liou, and R. G. Coleman, Comparative petrotectonic study of five Eurasian ultrahighpressure metamorphic complex, Int. Geol. Rev., 37, 191 - 211, 1995. Escartin, J., and G. Hirth, Nondilatant brittle deformation of serpentinites: Implications for Mohr-Coulomb theory and the strength of faults, J. Geophys. Res., 102, 2897 - 2913, 1997. Escartin, J., G. Hirth, and B. Evans, Strength of slightly serpentinized peridotites: Implications for the tectonics of oceanic lithosphere, Geology, 29, 1023 - 1026, 2001. Etheridge, M. A., V. I. Wall, and S. F. Cox, High fluid


X - 18

GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS
Kincaid, C., and I. S. Sacks, Thermal and dynamic evolution of the upper mantle in subduction zones, J. Geophys. Res., 102, 12,295 - 12,315, 1997. Kirby, S. H., S. Stein, E. A. Okal, and D. C. Rubie, Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere, Rev. Geophys., 34, 261 - 306, 1996. Krogh, E. G., C. W. Oh, and J. G. Liou, Polyphase and anticlockwise P-T evolution for the Franciscan eclogites and blueschists from Jenner, California, USA, J. Metamorph. Geol., 12, 121 - 134, 1994. Manning, C. E., Effect of sediments on aqueous silica transport in subduction zones, in Subduction: Top to Bottom, Geophys. Monogr. Ser., vol. 96, edited by G. E. Bebout et al., pp. 277 - 284, AGU, Washington, D.C., 1996. Martin, S., and P. Tartarotti, Polyphase HP metamorphism in the ophiolitic glaucophanites of the lower St. Marcel Valley (Aosta, Italy), Ofioliti, 15, 135 - 156, 1986. Maruyama, S., J. G. Liou, and M. Terabayashi, Blueschists and eclogites in the world and their exhumation, Int. Geol. Rev., 38, 485 - 594, 1996. Mattinson, J. M., Geochronology of high-pressure lowtemperature Franciscan metabasites: A new approach using the U-Pb system, in Blueschists and Eclogites, edited by B. W. Evans and E. H. Brown, Mem. Geol. Soc. Am., 164, 95 - 105, 1986. Mauler, A., G. Godard, and K. Kunze, Crystallographic fabrics of omphacite, rutile and quartz in Vende ?e ec logites ( Armorican Massif , F r ance). C onsequences for deformation mechanisms and regimes, Tectonophysics, 342, 81 - 112, 2001. McDowell, F. W., D. H. Lehmann, P. R. Gucwa, D. Fritz, and J. C. Maxwell, Glaucophane schists and ophiolites of the northern Coast Ranges: Isotopic ages and their tectonic implications, Geol. Soc. Am. Bull., 95, 1373 - 1382, 1984. Mei, S., and D. L. Kohlstedt, Influence of water on plastic deformation of olivine aggregates, 2, Dislocation c reep regime, J. Geophys. Res. , 10 5 , 21,471 - 21,481, 2000. Moore, D. E., Metamorphic history of a high-grade blueschist exotic block from the Franciscan Complex, California, J. Petrol., 25, 126 - 150, 1984. Moore, D. E., and M. C. Blake, New evidence for polyphase metamorphism of glaucophane schist and eclogite exotic blocks in the Franciscan Complex, California and Oregon, J. Metamorph. Geol., 7, 192 - 211, 1989. Moore, D. E., A. D. Lockner, M. Shengli, R. Summers, and J. D. Byerlee, Strengths of serpentinite gouges at elevated temperatures, J. Geophys. Res., 102, 14,787 - 14,801, 1997. Oh, C. W., and J. G. Liou, Metamorphic evolution of two different eclogites in the Franciscan Complex, Lithos, 25, 41 - 53, 1990. Parkinson, C. D., The origin and significance of metamorphosed tectonic blocks in melanges: Evidence from Sulsawesi, Indonesia, Terra Nova, 8, 312 - 323, 1996. Peacock, S. M., Creation and preservation of subduction related inverted metamorphic gradients, J. Geophys. Res., 92, 12,763 - 12,781, 1987a. Peacock, S. M., Serpentinization and infiltration metasomatism in the Trinity peridotite, Klamath provin c e, nort h ern C ali f ornia: Impli catio ns for subduction zones, Contrib. Mineral. Petrol., 95, 55 - 70, 1987b. Peacock, S., Numerical simulation of metamorphic pressure-temperature-time paths and fluid production in subduction slabs, Tectonics, 9, 1197 - 1211, 1990a. Peacock, S ., Fluid processes i n s ub du ct ion z on es, Science, 248, 329 - 337, 1990b. Peacock, S. M., The importance of blueschist ! eclogite dehydration reactions in the subduction oceanic crust, Geol. Soc. Am. Bull., 105, 684 - 694, 1993. Peacock, S. M., Thermal and petrologic structure of subduction zones, in Subduction: Top to Bottom, Geophys. Monogr. Ser., vol. 96, edited by G. E. Bebout et al., pp. 119 - 133, AGU, Washington, D.C., 1996. Perchuk, A. L., and P. Philippot, Geospeedometry and timescales of high-pressure m etamorphism, Int. Geol. Rev., 42, 207 - 223, 2000a. Perchuk, A. L., and P. Philippot, Nascent subduction: record in Yukon eclogites, Petrol ogy, 8, 1-16, 2000b. Perchuk, A. L., V. O. Yapaskurt, and S. K. Podlesskii, Genesis and exhumation dynamics of eclogites in the Kokchetav massif near mount Sulu-Tjube, Kazakhstan, Geochem. Int., 36, 877 - 885, 1998. Perchuk, A. L., P. Philippot, P. Erdmer, and M. Fialin, Rates of thermal equilibration at the onset of subduction deduced from diffusion modeling of eclogitic garnets, Yukon-Tanana terrain, Geology, 27, 531 - 534, 1999. Phillipot, P., Fluid-melt-rock interaction in mafic eclog ite s a nd c oe si t e b ea r i n g m e ta se d i m e nt s: Co nstraints on volatile recycling during subduction, Chem. Geol., 108, 93 - 112, 1993. Philippot, P., J. Blichert-Toft, A. L. Perchuk, S. Costa, and V. Yu. Gerasimov, Lu-Hf and Ar-Ar geochronology supports extreme rate of subduction zone metamorphism deduced from geospeedometry, Tectonophysics, 342, 23 - 38, 2001. Platt, J. P., Exhumation of high-pressure rocks: a review of concepts and processes, Terra Nova, 5, 119 - 133, 1993. Polino, R., G. V. Dal Piaz, and G. Gosso, Tectonic erosion at the Adria margin and accretionary processes for the Cretaceous orogeny of the Alps, ? ? Mem. Soc. Geol. Fr., 156, 345 - 367, 1990. Pysklywec, R. N., Evolution of subducting mantle lithosphere at a continental plate boundary, Geophys. Res. Lett., 28(23), 4399 - 4402, 2001. Pysklywec, R. N., C. Beaumont, and P. Fullsack, Modeling the behavior of the continental mantle lithosphere duri ng pl at e c onverg en ce, Geol ogy, 28, 655 - 658, 2000. Ranalli, G., Rheology of the Earth, 2nd ed., 413 pp., Chapman and Hall, New York, 1995. Regenauer-Lieb, K., D. A. Yuen, and J. Branlund, The initiation of subduction: Critically by addition of water?, Science, 294, 578 - 580, 2001. Reinecke, T., Very-high-pressure metamorphism and uplift of coesite-bearing sediments from the Zermatt-Saas zone, Western Alps, Eur. J. Mineral., 3, 7 - 17, 1991. Reinecke, T., Prograde high- to ultrahigh-pressure metamorphism and exhumation of oceanic sediments at Lago di Cignana, Zermatt-Saas Zone, western Alps, Lithos, 42, 147 - 189, 1998. Renner, J., B. Stockhert, A. Zerbian, K. Roller, and F. Å Rummel, An experimental study into the rheology of synthetic polycrystalline coesite aggregates, J. Geophys. Res., 106, 19,411 - 19,429, 2001. Ring, U., M. T. Brandon, S. D. Willett, and G. S. Lister, Exhumation processes, in Exhumation Processes: Normal Faulting, Ductile Flow, and Erosion, edited by U. Ring et al., Geol. Soc. Spec. Publ.,154, 1 - 27, 1999. Rubatto, D., and J. Hermann, Exhumation as fast as subduction, Geology, 29, 3 - 6, 2001. Ruff, L. J., and B. W. Tichelaar, What controls the seismogenic plate interface in subduction zones?, in Subduction: Top to Bottom, Geophys. Monogr. Ser., vol. 96, edited by G. E. Bebout et al., pp. 105 - 112, AGU, Washington, D.C., 1996. Scambelluri, M., and P. Philippot, Deep fluids in subduction zones, Lithos, 55, 213 - 227, 2001. Schmeling, H., R. Monz, and D. C. Rubie, The influence of olivine metastability on the dynamics of subduction, Earth Planet. Sci. Lett., 165, 55 - 66, 1999. Schmidt, M. W., and S. Poli, Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation, Earth Planet. Sci. Lett., 163, 361 - 379, 1998. Schott, B., and H. Schmeling, Delamination and detachment of a lithospheric root, Tectonophysics, 296, 225 - 247, 1998. Schreyer, W., Ultradeep metamorphic rocks: The retrospective viewpoint, J. Geophys., Res., 100, 8353 - 8366, 1995.

pressure during regional metamorphism and deformation: Implication for mass transport and deformation mechanism, J. Geophys. Res., 89, 4344 - 4358, 1984. Evans, B., and D. L. Kohlstedt, Rheology of rocks, in Rock Physics and Phase Relations, A Handbook of Physical Constants, AGU Ref. Shelf, vol. 3, edited by T. J. Ahrens, pp. 148 - 165, AGU, Washington, D.C., 1995. Gebauer, D., H. P. Schertl, M. Brix, and W. Schreyer, 35 Ma old ultrahigh-pressure metamorphism and evidence for very rapid exhumation in the Dora Maira Massif, Western Alps, Lithos, 41, 5 - 24, 1997. Gerya, T. V., and B. Stockhert, Exhumation rates of Å high pressure metamorphic rocks in subduct ion channels: The effect of rheology, Geophys. Res. Lett., 29(8), 1261, 10.1029/2001GL014307, 2002. Gerya, T. V., L. L. Perchuk, D. D. van Reenen, and C. A. Smit, Two-dimensional numerical modeling of pressure-temperature-time paths for the exhumation of some granulite facies terrains in the Precambrian, J. Geodyn., 30, 17 - 35, 2000. Grawinkel, A., and B. Stockhert, Hydrostatic pore fluid Å pressure to 9-km depth: Fluid inclusion evidence from the KTB deep drill hole, Geophys. Res. Lett., 24(24), 3273 - 3276, 1997. Green, H. W., and H. Houston, The mechanics of deep earthquakes, Annu. R ev. Earth Planet. Sci ., 23, 169 - 213, 1995. Guillot, S., K. H. Hattori, and J. de Sigoyer, Mantle wedge serpentinization and exhumation of eclogites: Insight from eastern Ladakh, northwest Himalaya, Geology, 28, 199 - 202, 2000. Guillot, S., K. H. Hattori, J. de Sigoyer, T. Nagler, and Å A. L. Auzende, Evidence of hydration of the mantle wedge and its role in the exhumation of eclogites, Earth Planet. Sci. Lett., 193, 115 - 127, 2001. Guinchi, C., and Y. Ricard, High-pressure/low-temperature metamorphism and the dynamic of an accretionary wedge, Geophys. J. Int., 136, 620 - 628, 1999. Hacker, B. R., Eclogite formation and the rheology, buoyancy, seismicity, and H2O content of oceanic crus t, in Su bdu ctio n: Top t o B ot t om , Ge op hy s. Monogr. Ser., vol. 96, edited by G. E. Bebout et al., pp. 337 - 346, AGU, Washington, D.C., 1996. Hacker, B. R., and S. M. Peacock, Creation, preservation, and exhumation of UHPM Rocks, in Ultrahigh Pressure Metamorphism, edited by R. G. Coleman and X . Wang, pp. 1 59 - 181, C ambr idge Un i v. Press, New York, 1995. Harley, S. L., and D. A. Carswell, Ultradeep crustal metamorphism: A prospective view, J. Geophys. Res., 100, 8367 - 8380, 1995. Hermann, J., O. Muntener, and M. Scambelluri, The importance of serpentinite mylonites for subduction and exhumation of oceanic crust, Tectonophysics, 327, 225 - 238, 2000. Hickman, S. H., and B. Evans, Kinetics of pressure solution at halite-silica interfaces and intergranular clay films, J. Geophys. Res., 100, 13,113 - 13,132, 1995. Hsu, K. J., Franciscan melange as a model for eugeosuncli nal sedimentation and underthrusting t ectonics, J. Geophys. Res., 76, 1162 - 1170, 1971. Huenges, E., J. Erzinger, J. Kuck, B. Engeser, and W. Kessels, The permeable crust: Geohydraulic properti es down to 9101 m depth, J. Geophys. Res., 102, 18,255 - 18,265, 1997. Karason, H., and R. D. van der Hilst, Tomographic imaging of the lowermost mantle with differential times of refracted and diffracted core phases (PKPPdiff), J. Geophys. Res., 106, 6569 - 6587, 2001. K as ah ar a, J . , A . K am i m ur a , G . F uj i e , an d R . H i no , Influence of water on earthquake generation along subduction zones, Bull. Earthquake Res. Inst. Univ. Tokyo, 76, 291 - 303, 2001. Kerrick, D. M., and J. A. D. Connolly, Metamorphic devolatilization of subducted oceanic metabasalts: Implications for seismicity, arc magmatism and volatile recycling, Earth Planet. Sci. Lett., 189, 19 - 29, 2001.


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS
Schwarz, S., and B. Stockhert, Pressure solution in siÅ liciclastic HP-LT metamorphic rocks--Constraints on the state of stress in deep levels of accretionary complexes, Tectonophysics, 255, 203 - 209, 1996. Schwartz, S., P. Allemand, and S. Guillot, Numerical model of the effect of serpentinization on the exhumation of eclogitic rocks: Insights from the Monviso ophilitic massif (Western Alps), Tectonophysics, 342, 193 - 206, 2001. Shen, A. H., and H. Keppler, Direct observation of complete miscibility the albite-H2O system, Nature, 385, 710 - 712, 1997. Shreve, R. L., and M. Cloos, Dynamics of sediment subduction, melange formation, and prism accretion, J. Geophys. Res., 91, 10,229 - 10,245, 1986. Sibson, R. H., Faulting and fluid flow, in Fluids in Tectonically Active Regimes of th e Continental Crust, edited by B. E. Nesbitt, pp. 92 - 132, Mineral. Assoc. of Can., Vancouver, 1990. Sobolev, N. V., and V. S. Shatsky, Diamond inclusions in garnets from metamorphic rocks: A new environment for diamond formation, Nature, 343, 742 - 746, 1990. Spakman, W., Delay-time tomography of the upper mantle below Europe, the Mediterranean, and Asia Minor, Geophys. J. Int., 107, 309 - 332, 1991. Stockhert, B., Stress and deformation in subduction Å zones - insight from the record of exhumed high pressure metamorphic rocks, in Deformation Mechanisms, Rheology, and Tectonics: Current Status and Future Perspectives, edited by S. de Meer et al., Geol. Soc. Spec. Publ., 200, 255 - 274, 2002. Stockhert, B., and J. Renner, Rheology of crustal rocks Å at ultrahigh-pressure, in When Continents Collide: Geodynamics and Geochemistry of Ultrahigh-Pressure Rocks, edited by B. R. Hacker and J. G. Liou, pp. 57 - 95, Kluwer Acad., Norwell, Mass., 1998. Suppe, J., and R. L. Armstrong, Potassium-argon dating of Franciscan metamorphic rocks, Am. J. Sci., 272, 217 - 233, 1972. Terabayashi, M., S. Maruyama, and J. G. Liou, Thermobaric structure of the Franciscan Complex in the Pacheco Pass Region, Diablo Range, California, J. Geol., 104, 617 - 636, 1996. Turcotte, D. L., and G. Schubert, Geodynamics: Application of Continuum Physics to Geological Problems, 450 pp., John Wiley, New York, 1982. van den Beukel, J., and R. Wortel, Thermo-mechanical model of arc-trench regions, Tectonophysics, 154, 177 - 193, 1988. van der Klauw, S. N. G. C., T. Reinecke, and B. Stockher t , E xh umat io n o f u lt rah i gh- pres su re met a m o rp hic o ce anic cr ust f r o m L ago d i C ig nan a , Piemontese zone, western Alps: The structural record in metabasites, Lithos, 41, 79 - 102, 1997. van Hunen, J., A. P. van den Berg, and N. J. Vlaar, A thermo-mechanical model of horizontal subduction below an overriding plate, Earth Planet. Sci. Lett., 182, 157 - 169, 2000. von Huene, R., and D. W. Scholl, Observations at con-

X - 19

vergent margins concerning sediment subduction, subduction erosion, and the growth of continental crust, Rev. Geophys., 29, 279 - 316, 1991. Wakabayashi, J., Counterclockwise P-T-t paths from amphibolites, Franciscan Complex, California: Relics from the early stages of subduction zone metamorphism, J. Geol., 98, 657 - 680, 1990. Weinberg, R. B., and H. Schmeling, Polydiapirs: multiwavelength gravity structures, J. Struct. Geol., 14, 425 - 436, 1992. Wheeler, J., Importance of pressure solution and coble creep in the deformation of polymineralic rock, J. Geophys. Res., 97, 4579 - 4586, 1992. Willett, S. D., Rheological dependence of extension in wedge models of convergent orogens, Tectonophysics, 305, 419 - 435, 1999. Yuan, X., et al., Subduction and collision processes in the Central Andes constrained by converted seismic phases, Nature, 408, 958 - 961, 2000.

À T. À Gerya and À V.À À À ÀÀ

B. Stockhert, Institut fur Geologie, Å Å Mineralogie und Geophysik, Ruhr-Universitat Bochum, Å D-44780 Bochum, Germany. (Taras.Gerya@ruhr-unibochum.de; Bernhard.Stoeckhert@ruhr-uni-bochum. de) A. L. Perchuk, Institute for Ore Deposits Geology, Petrology, Mineralogy, and Geochemistry, Russian Academy of Sciences, Staromonetny 35, 109017 Moscow, Russia. (alp@igem.ru)


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 4. Development of a wedge-shaped circulation pattern of material during ongoing subduction. The initial and boundary conditions correspond to 140 Â 100 km model A described in Table 1. Left: Evolution of the temperature field (isotherms labeled in œC) and distribution of rock types (color code): 1 = sedimentary rocks; 2 = basaltic upper oceanic crust; 3 = gabbroic lower oceanic crust; 4 = unhydrated mantle; 5 = serpentinized mantle; 6 = hydrated mantle beyond stability field of antigorite. Right: Evolution of viscosity (color code) and velocity field (arrows).

X-8


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 5. Development of a channel pattern of circulation of material during ongoing subduction. The initial and boundary conditions correspond to 140 Â 100 km Model B in Table 1. Colors and symbols correspond to those given in Figure 4.

X-9


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 6. Variations in geometry and thermal structure of the subduction zone predicted by the numerical models after 500 km of convergence, depending on the rate of subduction (a and b) and the characteristics of the hydration process (c and d), see text for details. The initial and boundary conditions correspond to models C (a), D (b), E (c) and F (d) described in Table 1. The inset separate sketches represent enlarged 110 Â 100 km areas of the original 140 Â 100 km models. Colors and symbols correspond to those given in Figure 4.

X - 10


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 7. Variations in the geometry and the thermal structure of the subduction zone after 500 km of convergence, depending on the initial age of the lithosphere (a and b) and on the inclination of the slab a (c and d). The initial and the boundary conditions correspond to models G (a), H (b), I (c) and J (d) specified in Table 1. The inset separate sketches represent enlarged 110/tg(a) Â 100 km areas of the original 140/tg(a) Â 100 km models. Colors and symbols correspond to those given in Figure 4.

X - 11


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 8. Variations in the geometry of the subduction zone after 1000 km of convergence, depending on the viscosity of serpentinized mantle. The initial and boundary conditions correspond to models L (a), and M (b) specified in Table 1. The separate sketches represent enlarged 110 Â 100 km areas of the original 140 Â 100 km models. Colors and symbols correspond to those given in Figure 4.

Figure 9. Particle trajectories and P-T paths calculated for the model with a wedge-shaped circulation pattern. The P-T paths represent the history of exhumed high-pressure metamorphic rocks. The initial and boundary conditions correspond to Model A specified in Table 1. The sketch in (a) represents the original 140 Â 100 km model geometry for the instant of ``sampling'' of the particle trajectories (b) and the corresponding P-T paths (c). Colors and symbols correspond to those given in Figure 4.

X - 12


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 10. Contrasting P-T evolution of two fragments of oceanic crust predicted for subduction with a moderate rate of 30 km/Myr. The initial and boundary conditions correspond to model A specified in Table 1. Separate sketches represent enlarged 74 Â 71 km areas of the original 140 Â 100 km model. Colors and symbols correspond to those given in Figure 4.

6 - 13


GERYA ET AL.: EXHUMATION OF HIGH-PRESSURE ROCKS

Figure 12. P-T evolution of two fragments of oceanic crust during subduction and exhumation. The initial and boundary conditions correspond to Model K specified in Table 1. The separate sketches represent enlarged 43 Â 71 areas of the original 81 Â 100 km model. Colors and symbols correspond to those given in Figures 4 and 9, respectively.

6 - 15