Numerical Geodynamic Modeling of Continental Convergent Margins

The continental convergence (subduction/collision) normally follows the oceanic subduction under the convergent forces of lateral ridge push and/or oceanic slab pull (Turcotte and Schubert, 2002). During these scenarios, a large amount of positively buoyant materials enter the trench causing slow down of the convergence that, eventually, may stop. However, before collision ceases, convergence between the plates can continue actively for tens of millions of years after ocean closure as it is testified by the 50 Ma active collisions in the Western Alps and Himalayas (e.g. Yin, 2006). A remarkable event during the early continental collision is the formation and exhumation of high-pressure to ultra-high-pressure (HP-UHP) metamorphic rocks, which is one of the most provocative findings in the Earth sciences during the past three decades. Occurrences of UHP terranes around the world have been increasingly recognized with more than 20 UHP terranes documented (e.g. Liou et al., 2004), which have repeatedly invigorated the concepts of deep subduction (>100 km) and subsequent exhumation of crustal materials (e.g. Chopin, 2003). It has been suggested that the HP-UHP metamorphism can be considered as a "hallmark" for the modern plate tectonics regime characterized by colder subduction and started from a Neoproterozoic time (e.g. Brown, 2006, 2007). The understanding of the dynamics of continental convergent margins implies several different but strictly correlated processes, such as continental deep subduction, HP-UHP metamorphism, exhumation, continental collision and mountain building. Besides the systematic geological/geophysical studies of the continental convergent zones, numerical modeling becomes a key and efficient tool (e.g. Burov et al., 2001; Yamato et al., 2007; Gerya et al., 2008; Warren et al., 2008a,b; Li and Gerya, 2009; Beaumont et al., 2009; Li et al., 2011). The tectonic styles of continental subduction can be either one-sided (overriding plate does not subduct) or two-sided (both plates subduct together) (Tao and O’Connell, 1992; Pope and Willett, 1998; Faccenda et al., 2008; Warren et al., 2008a), as well as several other possibilities, e.g. thickening, slab break-off, slab drips etc (e.g. Toussaint et al., 2004a,b). Models of HP-UHP rocks exhumation can be summarized into the following groups: (1) syn-collisional 13


Introduction
The continental convergence (subduction/collision) normally follows the oceanic subduction under the convergent forces of lateral ridge push and/or oceanic slab pull (Turcotte and Schubert, 2002).During these scenarios, a large amount of positively buoyant materials enter the trench causing slow down of the convergence that, eventually, may stop.However, before collision ceases, convergence between the plates can continue actively for tens of millions of years after ocean closure as it is testified by the 50 Ma active collisions in the Western Alps and Himalayas (e.g.Yin, 2006).A remarkable event during the early continental collision is the formation and exhumation of high-pressure to ultra-high-pressure (HP-UHP) metamorphic rocks, which is one of the most provocative findings in the Earth sciences during the past three decades.Occurrences of UHP terranes around the world have been increasingly recognized with more than 20 UHP terranes documented (e.g.Liou et al., 2004), which have repeatedly invigorated the concepts of deep subduction (>100 km) and subsequent exhumation of crustal materials (e.g.Chopin, 2003).It has been suggested that the HP-UHP metamorphism can be considered as a "hallmark" for the modern plate tectonics regime characterized by colder subduction and started from a Neoproterozoic time (e.g.Brown, 2006Brown, , 2007)).The understanding of the dynamics of continental convergent margins implies several different but strictly correlated processes, such as continental deep subduction, HP-UHP metamorphism, exhumation, continental collision and mountain building.Besides the systematic geological/geophysical studies of the continental convergent zones, numerical modeling becomes a key and efficient tool (e.g.Burov et al., 2001;Yamato et al., 2007;Gerya et al., 2008;Warren et al., 2008a,b;Li and Gerya, 2009;Beaumont et al., 2009;Li et al., 2011).The tectonic styles of continental subduction can be either one-sided (overriding plate does not subduct) or two-sided (both plates subduct together) (Tao and O'Connell, 1992;Pope and Willett, 1998;Faccenda et al., 2008;Warren et al., 2008a), as well as several other possibilities, e.g.thickening, slab break-off, slab drips etc (e.g.Toussaint et al., 2004a,b).Models of HP-UHP rocks exhumation can be summarized into the following groups: (1) syn-collisional exhumation of a coherent and buoyant crustal slab, with formation of a weak zone at the entrance of the subduction channel (Chemenda et al., 1995(Chemenda et al., , 1996;;Toussaint et al., 2004b;Li and Gerya, 2009); (2) episodic ductile extrusion of HP-UHP rocks from the subduction channel to the surface or crustal depths (Beaumont et al., 2001;Warren et al., 2008a); (3) continuous material circulation in the rheologically weak subduction channel stabilized at the plate interface, with materials exhumed from different depths (Burov et al., 2001;Stöckhert and Gerya, 2005;Yamato et al., 2007;Warren et al., 2008a).In this chapter, the processes and dynamics of continental subduction/collision and HP-UHP rocks exhumation are investigated by the method of large-scale numerical geodynamic modeling.First the numerical method is described, which is followed by the numerical model setup and systematic thermo-mechanical numerical experiments.The discussion section covers a broad range of topics related to the continental subduction and exhumation.Finally a concluding part is presented.

Governing equations and numerical implementation
The momentum, continuity and heat conservation equations for a 2D creeping flow including thermal and chemical buoyant forces are solved: where the density ρ depends on composition (C), melt fraction (M), pressure (P) and temperature (T); g is the acceleration due to gravity.
(ii) Conservation of mass is approximated by the incompressible continuity equation q x = −k(C, P, T) ∂T ∂x , q z = −k(C, P, T) ∂T ∂z xz εxz where D/Dt is the substantive time derivative.x and z denote the horizontal and vertical directions, respectively.The deviatoric stress tensor is defined by σ ′ xx , σ ′ xz , σ ′ zz , whilst the strain rate tensor is defined by εxx , εxz , εzz .q x and q z are heat flux components.ρ is the density.k(C, P, T) is the thermal conductivity as a function of composition (C), pressure (P) and temperature (T).C p is the isobaric heat capacity.H r , H a , H s , H L are radioactive, adiabatic, shear and latent heat production, respectively (see Table 1 for details of these parameters).To solve the above equations, the I2VIS code is used (Gerya and Yuen, 2003a).It is a two-dimensional finite difference code with marker-in-cell technique which allows for non-diffusive numerical simulation of multi-phase flow in a rectangular fully staggered Eulerian grid.I2VIS accounts for visco-plastic deformation and several geological processes that are described below.All abbreviations and units used in this chapter are listed in

Boundary conditions
For the 2D numerical models presented in this chapter, the velocity boundary conditions are free slip at all boundaries except the lower one, which is permeable (Burg and Gerya, 2005; 275 Numerical Geodynamic Modeling of Continental Convergent Margins www.intechopen.comLi et al., 2010).Infinity-like external free slip conditions along the lower boundary imply free slip to be satisfied at 1000 km below the bottom of the model.As for the usual free slip condition, external free slip allows global conservation of mass in the computational domain and is implemented by using the following limitation for velocity components at the lower boundary: ∂v x /∂z = 0, ∂v z /∂z = −v z /∆z external , where ∆z external is the vertical distance from the lower boundary to the external boundary where free slip (∂v x /∂z = 0, v z = 0) is satisfied.
The thermal boundary conditions have a fixed value (0 • C) for the upper boundary and zero horizontal heat flux across the vertical boundaries.For the lower thermal boundary, an infinity-like external constant temperature condition is imposed, which allows both temperatures and vertical heat fluxes to vary along the permeable box lower boundary, implying constant temperature condition to be satisfied at the external boundary.This condition is implemented by using the limitation ∂T/∂z =( T external − T)/∆z external where T external is the temperature at the external boundary and ∆z external is the vertical distance from the lower boundary to the external boundary (Burg and Gerya, 2005;Li et al., 2010).

Rheological model
A viscoplastic rheology is assigned for the model in which the rheological behaviour depends on the minimum differential stress attained between the ductile and brittle fields.
Ductile viscosity dependent on strain rate, pressure and temperature is defined in terms of deformation invariants as: where ε II = 0.5 εij εij is the second invariant of the strain rate tensor.A D , E, V and n are experimentally determined flow law parameters (Table 2).F is a dimensionless coefficient depending on the type of experiments on which the flow law is based.For example: F = [2 (1−n)/n ]/[3 (1+n)/2n ] for triaxial compression and F = 2 (1−2n)/n for simple shear.
The ductile rheology is combined with a brittle/plastic rheology to yield an effective visco-plastic rheology.For this purpose the Mohr-Coulomb yield criterion (e.g.Ranalli, 1995) is implemented as follows: 2 ε II where σ yield is the yield stress.ε II is the second invariant of the strain rate tensor.P is the dynamic pressure.C is the cohesion.ϕ is the internal frictional angle.λ is the pore fluid coefficient that controls the brittle strength of fluid-containing porous or fractured media (Brace and Kohlstedt, 1980).ϕ ef f can be illustrated as the effective internal frictional angle that integrates the effects of internal frictional angle (ϕ) and pore fluid coefficient (λ).λ is the pore fluid coefficient that controls the brittle strength of fluid-containing porous or fractured media.The effective viscosity of molten rocks (M ≥ 0.1) was calculated using the formula (Pinkerton and Stevenson, 1992;Bittner and Schmeling, 1995): where η 0 is an empirical parameter depending on rock composition, being η 0 = 10 13 Pa s (i.e. 1 × 10 14 ≤ η ≤ 2 × 10 15 Pa s for 0.1 ≤ M ≤ 1) for molten mafic rocks and η 0 = 5 × 10 14 Pa s (i.e. 6 × 10 15 ≤ η ≤ 8 × 10 16 Pa s for 0.1 ≤ M ≤ 1) for molten felsic rocks.Successfully tested for a broad range of suspensions with various bubble or crystal conventions, this formula takes into account, other than concentration, particle shape and size distribution.

Partial melting model
The numerical code accounts for partial melting of the various lithologies by using experimentally obtained P-T dependent wet solidus and dry liquidus curves (Gerya and Yuen, 2003b).As a first approximation, volumetric melt fraction M is assumed to increase linearly with temperature accordingly to the following relations (Burg and Gerya, 2005): where T solidus and T liquidus are the wet solidus and dry liquidus temperature of the given lithology, respectively (Table 3).Consequently, the effective density, ρ ef f , of partially molten rocks varies with the amount of melt fraction and P-T conditions according to the relations: where ρ solid and ρ molten are the densities of the solid and molten rock, respectively, which vary with pressure and temperature according to the relation: where ρ 0 is the standard density at P 0 = 0.1 MPa and T 0 = 298 K; α and β are the thermal expansion and compressibility coefficients, respectively (Tables 1 and 3).The effects of latent heat H L (e.g.Stüwe, 1995) are accounted for by an increased effective heat capacity (C Pe f f ) and thermal expansion (α ef f ) of the partially molten rocks (0 < M < 1), calculated as where C P and α are the heat capacity and the thermal expansion of the solid crust, respectively, and Q L is the latent heat of melting of the crust (Table 1).

Topographic model
The spontaneous deformation of the upper surface of the lithosphere, i.e. topography, is calculated dynamically as an internal free surface by using a low viscosity (e.g. 10 18 Pa s), initially 8-12 km thick layer (thickness of this layer changes dynamically during experiments) above the upper crust.The composition is either "air" (1 kg/m 3 , above water level) or "water" (1000 kg/m 3 , below water level).The interface between this weak layer and the underlying crust is treated as an internal erosion/sedimentation surface which evolves according to the Eulerian transport equation solved in Eulerian coordinates at each time step (Gerya and Yuen, 2003b): where z es is the vertical position of the surface as a function of the horizontal distance v x .v z and v x are the vertical and horizontal components of the material velocity vector at the surface.v s and v e are the sedimentation and erosion rates, respectively, which correspond to the relation: v s = 0, v e = v e0 , when z es < erosion level; v s = v s0 , v e = 0, when z es > erosion level; where v e0 and v s0 are the imposed constant large scale erosion and sedimentation rates, respectively.The code allows for marker transmutation that simulates erosion (rock markers are transformed to weak layer markers) and sedimentation (weak layer markers are transformed to sediments).
The partially molten crustal rocks (13,14,15,16,17,18) are not shown in this figure, which will appear during the evolution of the model.In our numerical models, the medium scale layering usually shares the same physical properties, with different colors used only for visualizing plate deformation.Detailed properties of different rock types are shown in Tables 2 and 3.

279
Numerical Geodynamic Modeling of Continental Convergent Margins www.intechopen.comconvergence velocity (V x ) in a small internal domain that remains fixed with respect to the Eulerian coordinate (Fig. 1).
In the numerical models, the driving mechanism of subduction is a combination of plate push (prescribed rightward convergence velocity) and increasing slab pull (temperature-induced density contrast between the subducted lithosphere and surrounding mantle).This type of boundary condition is commonly used in numerical models of subduction and collision (e.g.Toussaint et al., 2004b;Burg and Gerya, 2005;Currie et al., 2007;Yamato et al., 2007;Warren et al., 2008b) and assumes that in the globally confined three-dimensional system of plates, local external forcing coming either from different slabs or from different sections of the same laterally non-uniform slab can be significant.
Following previous numerical studies performed with similar geodynamic settings (e.g. Warren et al., 2008a;Li and Gerya, 2009) we design numerical models consisting of three major domains (from left to right, Fig. 1): (1) a pro-continental domain, (2) an oceanic domain, and (3) a retro-continental domain.The subducting pro-continent comprises a marginal unit and an interior unit.In the continental domain, the initial material field is set up by a 35 km thick continental crust composed of sediment (6 km thick), upper crust (14 km thick) and lower crust (15 km thick), overlying the lithospheric mantle (85 km thick) and subjacent mantle (540 km thick).The oceanic domain comprises an 8 km thick oceanic crust overlying the lithospheric mantle (82 km thick) and subjacent mantle (570 km thick).The material properties of all layers (Fig. 1) are listed in Table 3.The initial thermal structure of the lithosphere (white lines in Fig. 1) is laterally uniform with 0 • C at the surface and 1300 • C at the bottom of the lithospheric mantle (both continental and oceanic).The initial temperature gradient in the asthenospheric mantle is around 0.5 • C/km.

Model result 4.1 Reference model
The reference model is designed with prescribed convergence velocity (V x )of5cm/y.All the other configurations and parameters are shown in Figure 1 and Table 3.
At the initial stages, the relatively strong oceanic plate subducts along the weak zone to the mantle (Fig. 2a).The continental margin subducts to >100 km depth, following the high-angle oceanic subduction channel (Fig. 2b).The significant characters are the detachment of subducting upper/middle crust at the entrance zone of the subduction channel with a series of thrust faults formed (Fig. 2b-d).A small amount of crustal rocks located in the lower part of the channel are detached from the plate at asthenospheric depths, indenting into the mantle wedge and forming a compositionally buoyant plume (Fig. 2c).Such sub-lithospheric plumes are discussed in detail in Currie et al. ( 2007) and Li and Gerya (2009).In addition, a partially molten plume forms in the deeply subducted oceanic plate and moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 2c,d).As subduction continues, another partially molten plume forms in the deeply subducted continental plate.It also moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 2e,f).The characteristics and 2D and 3D dynamics of this kind of plume are studied in detail in Gerya and Yuen (2003b) and Zhu et al. (2009).
As continental subduction continues, partially molten rocks accumulated in the subduction channel extrude upward to the crustal depths (Fig. 2d,e).Then these UHP rocks exhume buoyantly to the surface forming a dome structure (Fig. 2f).The exhumed UHP rocks are mainly located near the suture zone with a fold-thrust belt formed in the foreland extending for about 300-400 km (Fig. 2f).P-T paths (Fig. 2, inset) show that peak P-T conditions 280 Earth Sciences www.intechopen.com of the exhumed rocks are 2.5-4 GPa and 600-800 • C, respectively (also see Figure 3 for the peak pressure and temperature conditions of the collision zone).This indicates the UHP metamorphic rocks are formed and exhumed from a depth >100 km.

Models with variable convergence velocity
The reference model is further investigated with lower convergence velocity (2.5 cm/y) and higher convergence velocity (10 cm/y).All the other parameters are the same as in Tables 2  and 3.In the slow convergence regime, the continental margin also subducts to >100 km depth along the high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 4a).
The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 4a-d).With continued continental subduction, partially molten rocks accumulated in the subduction channel extrude upward to the crustal depth (Fig. 4c,d).P-T paths (Fig. 4) show that peak P-T conditions of the exhumed rocks are 2-4 GPa and 600-800 • C. In the fast convergence regime, the continental domain continues subducting along the high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 5a).
Then the crustal rocks in the lower part of the channel detach from the plate at asthenospheric depths, intrude into the mantle wedge, and form a horizontal compositionally buoyant plume (Fig. 5a-c).In addition, a partially molten plume forms in the deeply subducted plate and moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 5c,d), which is similar to behavior of the reference model.After the convergence ceases (1500 km shorting, 15 Myr), the subducted continental crustal rocks in the sub-lithospheric plume extrude upward to the surface forming a dome structure (Fig. 5d).P-T paths (Fig. 5) indicate that peak P-T conditions of the exhumed rocks are 3-4 GPa and 600-800 • C, respectively.This parameter sensitivity studies indicate that the slower convergence produces very small sub-lithospheric plume (Fig. 4a), coupled subduction channel and wide collision zone (Fig. 4d).In contrast, the faster convergence results in very large sub-lithospheric plume (Fig. 5a), decoupled subduction channel and narrow collision zone (Fig. 5d).Both of the models can obtain UHP rocks exhumation.However, the convergence velocity changes the amount of crustal rocks subducted to and exhumed from UHP depth (c.f.Figs 4 and 5).The lithospheric thermal structure plays an important role on subduction/collision processes (e.g.Toussaint et al., 2004a, b).Therefore we investigate the sensitivity of oceanic thermal gradient for the reference model.In the hot model, initial thermal structure of the oceanic lithosphere is linearly interpolated with 0 • C at the surface (≤ 10 km depth) and 1300 • Cat70 km depth (compared to 1300 • C at 100 km depth in the reference model).In contrast, the initial thermal structure of oceanic lithosphere in the cold model is linearly interpolated with 0 • Cat the surface (≤ 10 km depth) and 1300 • Cat130km depth.In the hot model, the continental margin subducts following the oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 6a).Then the crustal rocks in the lower part of the channel detach from the plate at asthenospheric depths, intrude into the mantle wedge, and form a horizontal compositionally buoyant plume (Fig. 6b).The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 6b,c).With continued continental subduction, partially molten rocks in the middle channel extrude upward to the crustal depth (Fig. 6c,d).The subduction channel is highly coupled.As a result, the partially molten rocks in the sub-lithospheric plume stay at the bottom of the overriding lithosphere (without exhumation).In the cold model, the continental margin subducts following the oceanic plate to the bottom of the lithospheric mantle (Fig. 7a).In this case, there is no sub-lithospheric plume formed.Consequently, the subduction channel is thicker and thicker.The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 7b,c).The subducted continental crustal rocks extrude upward to the surface forming a dome structure (Fig. 7c,d).The subduction channel is highly decoupled.The shape and characteristics of the subduction channel in the hot model (Fig. 6) is similar to that in the slow convergence model (Fig. 4).It indicates that both the higher temperature and the slower convergence can increase the rheological coupling at plate interface.As a result, coupled subduction channel is produced in these two models.In contrast, decoupled channels are formed in the colder model (Fig. 7) as well as in the faster convergence model (Fig. 5).

Flow modes in the subduction channel
To a first approximation, viscous channel flow can be analysed using lubrication theory (e.g.England and Holland, 1979;Cloos, 1982;Cloos andShreve, 1988a, 1988b;Mancktelow, 1995; 285 Numerical Geodynamic Modeling of Continental Convergent Margins www.intechopen.comRaimbourg et al., 2007;Warren et al., 2008b;Beaumont et al., 2009).Under the lubrication approximations, channel flow velocity is where η is the assumed uniform viscosity in the subduction channel, ∂P/∂x is the effective down-channel pressure gradient, with x measured in the down-dip direction, y is the position in the channel measured normal to the base, h is channel thickness and U is the subduction velocity of the underlying lithosphere (Fig. 8a).The overlying lithosphere is assumed to be stationary.When nondimensional variables u ′ = u/U, h ′ = h/H, x ′ = x/h and y ′ = y/h are used, Equation 13 reduces to where η ef f , the effective viscosity and ∂P/∂x are averaged quantities measured at the channelscale used to estimate H.This parameter H is the characteristic channel thickness for E ∼ 1, the balance point between downward and return flows.Numerical models of the subduction channels can be conveniently interpreted in terms of the characteristic exhumation number E (Raimbourg et al., 2007;Warren et al., 2008b;Beaumont et al., 2009).The corresponding flow modes associated with burial and exhumation of UHP rocks are shown in Figure 8.The first order dynamics can be approximated by the above-mentioned lubrication theory for creeping flows, and characterized in terms of the competition between down-channel Couette flow (U(1 − y/h) in Eq. 13), caused by the drag of the subducting lithosphere and the opposing up-channel Poiseuille flow ((1/2 η)(∂P/∂x)(yh− y 2 ) in Eq. 13), driven by the buoyancy of low-density subducted crustal material.This competition can be expressed through the exhumation number E, which is a force ratio derived from the non-dimensional channel flow equation (Eq.14).The actual values that determine E (Eqs 15, 16) will depend on the particular problem and its evolving solution.For a channel with deformable walls and no tectonic over/under-pressure, ∂P/∂x =(ρ m − ρ c ) g sin(θ) where θ is channel dip (Fig. 8a).
Along with the characteristic E, defined at the scale of the subduction channel, space-time variations in the channel flow can be interpreted in terms of the local exhumation number E(x, t) and corresponding flow modes (Warren et al., 2008b;Beaumont et al., 2009).During continental subduction, E(x, t) evolves from <1 during subduction (c.f.Fig. 8b and Fig. 2a), to ∼1 during detachment and stagnation in the subduction channel (c.f.Fig. 8c and Fig. 2b,c), to >1 at the onset of and during exhumation (c.f.Fig. 8d and Fig. 2d-e).Buoyancy is a necessary, but not sufficient, condition for UHP exhumation.Among other controlling factors (Fig. 8), decreasing viscosity (η ef f ) is typically most important for driving E beyond the exhumation threshold.In general, E(x, t) should be regarded as a measure of local exhumation potential, even where the local threshold value is exceeded (E>1), efficient exhumation may be impeded by constrictions (small h) or high viscosities (large η ef f ) further up the channel.

Coupled and decoupled subduction channel
The numerical results show that the coupled subduction channel favors lower convergence velocity (Fig. 4) and hotter oceanic lithosphere (Fig. 6).It is characterized by continuous accretion of the weak upper continental crust resulting in the development of a thick and broad crustal wedge.In contrast, the higher convergence velocity (Fig. 5) and colder oceanic lithosphere (Fig. 7) result in decoupling of the convergent plates.Transition from coupled to decoupled regime occurs always at the early stages of continental collision indicating that insertion of rheologically weak crustal material in the subduction channel is critical for the subsequent evolution of the collision zone (Faccenda et al., 2009).The numerical models confirm that HP-UHP complexes can be formed in both coupled and decoupled channels in the wide range of convergence scenarios (Figs 2-7).
As discussed in Faccenda et al. (2009), coupled collision zones (which can be either retreating or advancing) are characterized by a thick crustal wedge and compressive stresses (i.e.Himalaya and Western Alps), while decoupled end-members (which are always retreating)

287
Numerical Geodynamic Modeling of Continental Convergent Margins www.intechopen.comare defined by a thin crustal wedge and bi-modal distribution of stresses (i.e.compressional in the foreland and extensional in the inner part of the orogen, Northern Apennines).

Thrust fault formation and exhumation of (U)HP units
Fig. 9. Highly-compressional regime of continental subduction (low pull force) after Chemenda et al. (1995Chemenda et al. ( , 2000)).One of the most important characteristics of the numerical models presented in this study is the formation of the thrust faults (rheological weak zones), which is followed by exhumation processes (e.g.Fig. 2).Similar detachment phenomenon is also documented in analogue models of continental subduction (Fig. 9; Chemenda et al., 1995Chemenda et al., , 1996Chemenda et al., , 2000)).
The behavior of the subducted continental crust depends on two competing effects: upward buoyancy and downward subduction drag as discussed in Section 5.1.Subduction drag within the crust and mantle drives the subduction of buoyant crustal materials into larger depths (Fig 2a).At the same time the buoyancy forces and also the deviatoric stresses increase in the subduction channel.When the materials are no longer strong enough to sustain the accumulated buoyancy and deviatoric stresses, the subducted continental crust will yield with forming the rheological weak zone (thrust fault) (Fig. 2b,c) followed by the detachment and exhumation of the buoyant crustal materials (Fig. 2d-f) and release of the accumulated buoyancy and deviatoric stresses.

Upper crustal structure of the HP-UHP terrane
The upper-crustal settings of many UHP terranes share a number of structural characteristics (Fig. 10a; Beaumont et al., 2009): (1) a dome structure cored by the UHP nappe, (2) domes flanked by low-grade accretionary wedge and/or upper crustal sedimentary rocks, (3) overlying and underlying medium-to high-pressure nappes, (4) suture zone ophiolites and (5) foreland-directed thrust faults and the syn-exhumation normal faults.Our numerical models reproduce the general characteristic upper crustal structures (Fig. 10b), especially the dome structure of the HP-UHP cores, the flanked and overlaid low-grade accretionary wedge,  the foreland-directed thrust faults and the fold-thrust belt.The ophiolites are distributed near the suture zone in the reference model (Fig. 10b).

Influence of tectonic overpressure on the P-T paths of (U)HP rocks
The principle of lithostatic pressure is habitually used in metamorphic geology to calculate burial/exhumation depth from pressure given by geobarometry.However pressure deviation from lithostatic, i.e. tectonic overpressure/underpressure due to deviatoric stress and deformation, is an intrinsic property of flow and fracture in all materials, including rocks under geological conditions (e.g.Rutland, 1965;Brace et al., 1970;Mancktelow, 1993Mancktelow, , 1995Mancktelow, , 2008;;Petrini and Podladchikov, 2000).Therefore, one important question is whether the principle of lithostatic pressure is applicable in subduction/collision zones where crystallization and exhumation of HP-UHP rocks take place.Some authors have argued that rocks under geological conditions are too weak to support significant overpressure (e.g.Brace et al., 1970;Ernst, 1971;Burov et al., 2001;Renner et al., 2001;Green, 2005).Yet, St üwe and Sandiford (1994) suggested that petrologically derived P-T paths may not record depth changes only but stress changes.In addition, lithospheric-scale numerical models reveal regions where pressure may be hundreds of MPa or even several GPa higher or lower than lithostatic values (Mancktelow, 1993(Mancktelow, , 1995;;Petrini and Podladchikov, 2000;Burg and Gerya, 2005).The analytical solutions show that the tectonic overpressure can be as high as 60% of the lithostatic value in the brittle regime.In contrast, the ductile overpressure is normally <10% of the lithostatic value (Li et al. 2010).Li et al. (2010) conducted systematic numerical simulations of continental subduction/collision zones with variable brittle and ductile rheologies of the crust and mantle.In the numerical model (Figs 11,12), the uppermost lithospheric mantle that can be considered as the wall of the subduction channel shows the largest tectonic overpressure (≥1 GPa and ≥50% of the lithostatic pressure).However, these overpressured zones rarely  Colours of these squares are used for discrimination of marker points plotted in the P-T diagrams and do not correspond to the colours of rock types.The solid curves show P-T paths with dynamic pressure, while the dashed curves show that with lithostatic pressure.
or never participate in the exhumation processes.Hence, they do not influence the P-T conditions of geologically distributed HP-UHP rocks in nature.The main overpressure region that may influence the P-T paths of HP-UHP rocks is located in the bottom corner of the wedge-like confined channel (Fig. 11) with the characteristic magnitude of pressure deviation on the order of ∼0.3 GPa and 10-20% from the lithostatic values (Fig. 12).The degree of confinement of the subduction channel is the key factor controlling this magnitude.The models show that the overpressure is small (∼10% lithostatic) and should not affect in a crucial way the metamorphic mineral equilibria of the exhumed UHP rocks.The challenge would be to identify the geological record to actually measure precisely such minor deviations.An important unresolved issue concerns the question of how the tectonic overpressure could be potentially recorded by mineral equilibria in natural rocks.Thermodynamic tools used for thermobarometry of natural rocks are mainly based on experimental data obtained under conditions of an isotropic stress state (i.e. in the absence of significant deviatoric stresses) and may not be directly applicable for recording dynamic pressure in strongly stressed rocks.Indeed, not all overpressured rocks should be strongly internally stressed.For example, weak (e.g.reacting, fluid rich) rock inclusions in strong overpressured stressed lithosphere will also have strong overpressure, but the stress state will be isotropic, i.e. similar to the conditions of laboratory experiments.This isotropic stress (dynamic pressure) will be notably different from the lithostatic pressure which will be then directly recorded by mineral equilibria of such rock inclusions.Obviously further efforts are needed to experimentally study mineral equilibria in stressed rocks.The tectonic overpressure is also investigated for several different subduction/collision and exhumation scenarios (Fig. 13).In these numerical models, the overpressures are not significant in the mature subduction channel and/or the inner collision belt, which suggests that the overpressure that can possibly affect the HP-UHP rocks is mainly related to the wedge-like confined subduction channel (Figs 11,12).

Conclusions and future perspectives
Continental collision was investigated with numerical models where an advancing oceanic/continental plate subducts under a fixed continent.The most important results can be summarized as follows: (1) During the processes from continental subduction to exhumation, the flow mode in the subduction channel changes from Couette-dominated to Poiseuille-dominated flows.Numerical models of the subduction channels can be conveniently interpreted in terms of the characteristic exhumation number.
(2) The coupled subduction channel is formed in the models with slower convergence or hotter oceanic lithosphere.Whereas faster convergence or colder oceanic lithosphere result in decoupling of the converging plates.
(3) The continental margin can subduct following the oceanic plate to >100 km depth, and then exhume to the surface forming a HP-UHP dome.The upper crustal structures of the collision zone in the numerical models are consistent with both the analogue model and the natural UHP zones.The exhumation of UHP rocks occurs in a large variety of numerical parameters.(4) The main tectonic overpressure region that may influence the P-T paths of HP-UHP rocks is located in the bottom corner of the wedge-like confined channel with the characteristic magnitude of pressure deviation on the order of ∼ 0.3 GPa and 10-20% from the lithostatic values.The degree of confinement of the subduction channel is the key factor controlling this magnitude.The tectonic overpressure should not affect in a crucial way the metamorphic mineral equilibria of the exhumed UHP rocks.The challenge would be to identify the geological record to actually measure precisely such minor deviations.The numerical modeling method is demonstrated to be a great tool to study the geodynamic processes in the continental convergent margins.Most of the existed numerical models in the relevant topics are based on the two-dimensional regimes with the along-strike variations ignored.However, the tectonic processes remain inherently three-dimensional.So it is quite significant to conduct 3D numerical modeling to investigate the dynamics of continental collision with applicable to the natural tectonic settings (e.g.Alpine and Himalayan collision belts).

Fig. 2 .Fig. 3 .
Fig. 2. Enlarged domain evolution (1300 × 600 km) of the reference model.Colors of rock types are as in Figure 1.Time (Myr) of shortening is given in the figures.White numbered lines are isotherms in • C. Small colored squares indicate positions of representative markers (rock units) for which P-T paths are shown (inset).Colors of these squares are used for discrimination of marker points plotted in P-T diagrams and do not correspond to the colors of rock types.

Fig. 4 .
Fig. 4. Enlarged domain evolution (800 × 225 km) of the model with lower convergence velocity V x = 2.5 cm/y.Colors of rock types are as in Figure 1.Time (Myr) of shortening is given in the figures.White numbered lines are isotherms in • C. Small colored squares indicate positions of representative markers (rock units) for which P-T paths are shown (right).Colors of these squares are used for discrimination of marker points plotted in P-T diagrams and do not correspond to the colors of rock types.

Fig. 5 .
Fig. 5. Enlarged domain evolution (1075 × 275 km) of the model with higher convergence velocity V x = 10 cm/y.Colors of rock types are as in Figure 1.Time (Myr) of shortening is given in the figures.White numbered lines are isotherms in • C. Small colored squares indicate positions of representative markers (rock units) for which P-T paths are shown (right).Colors of these squares are used for discrimination of marker points plotted in P-T diagrams and do not correspond to the colors of rock types.

Fig. 6 .
Fig. 6.Enlarged domain evolution (900 × 225 km) of the model with higher temperature of the oceanic lithosphere (hot model).Colors of rock types are as in Figure 1.Time (Myr)of shortening is given in the figures.White numbered lines are isotherms in • C. Small colored squares indicate positions of representative markers (rock units) for which P-T paths are shown (right).Colors of these squares are used for discrimination of marker points plotted in P-T diagrams and do not correspond to the colors of rock types.

Fig. 7 .
Fig. 7. Enlarged domain evolution (900 × 225 km) of the model with lower temperature of the oceanic lithosphere (cold model).Colors of rock types are as in Figure 1.Time (Myr)of shortening is given in the figures.White numbered lines are isotherms in • C. Small colored squares indicate positions of representative markers (rock units) for which P-T paths are shown (right).Colors of these squares are used for discrimination of marker points plotted in P-T diagrams and do not correspond to the colors of rock types.

Fig. 8 .
Fig. 8. Schematic diagram showing subduction/exhumation channel flow behavior in terms of dominating Couette (subduction) and Poiseuille (exhumation) flows (after Warren et al., 2008b; Beaumont et al., 2009).(a) General nomenclature.(b-d) Flow types identified in the models.(b) Couette flow (subduction) dominates.All flow is directed downward.(c) Buoyant materials stagnate at bottom of channel with the Poiseuille flow effect increasing.It is characterized as the transition from subduction-dominated to exhumation-dominated channel (d) Poiseuille flow (exhumation) dominates.Buoyancy-driven exhumation starts at channel bottom and propagates upward.

Fig. 11 .
Fig. 11.Evolution of the wedge-like subduction channel within enlarged 530 × 210 km domain of the original 4000 × 670 km model.Time (Myr) of shortening is given in the figures.White numbered lines are isotherms in • C. Small coloured squares (with '+' in them) show positions of representative markers (rock units) for which P-T paths are shown (right).Colours of these squares are used for discrimination of marker points plotted in the P-T diagrams and do not correspond to the colours of rock types.The solid curves show P-T paths with dynamic pressure, while the dashed curves show that with lithostatic pressure.

Table 1 .
Abbreviations and units of the variables used in this chapter.