## Abstract

High-pressure–low-temperature (HP-LT) metamorphic rocks that belong to the same orogen commonly show alignment of their peak pressure and related temperature within a pressure-temperature (PT) diagram, defining a pressure/temperature (P/T) ratio. In the Aegean region, for example, two metamorphic belts of different ages, the Eocene Cycladic blueschists and the Miocene Cretan blueschists, show contrasting P/T ratios that are characterized by the slope and the temperature at 10 kbar (ΔP/ΔT and T@10 kb). Peak PT data for the Cyclades yield a ΔP/ΔT of 0.056 and T@10 kb of 375 °C, while the Cretan blueschists give a much larger ΔP/ΔT with a lower T@10 kb. One-dimensional (1D) modeling of the thermal evolution of a subducted continental margin shows that subduction velocity controls both the slope ΔP/ΔT and T@10 kb of the P/T ratio, and the subduction dip only controls T@10 kb. On these bases, the variations of P/T ratios in the Aegean region reflect variations through time of subduction dip and velocity. Eocene subduction for Cycladic blueschists burial occurred at a rate of 1.5 cm a^{−1}, while subduction velocity during Cretan blueschists formation is found to be 2.75 cm a^{−1}. Because the convergence rate between Africa and Eurasia is constant and ~1–1.5 cm a^{−1} at these times, the active southward rollback of the Aegean slab during the Oligo-Miocene likely explains the larger subduction velocity for the Cretan HP-LT rocks. These results exemplify the use of this new modeling approach as a proxy to quantify dip and velocity of continental subduction from the P/T ratio of high-pressure rocks.

## INTRODUCTION

High-pressure–low-temperature (HP-LT) metamorphic rocks that belong to the same orogen commonly show alignment of their peak pressure and related temperature within a PT diagram. Such linear relationship between peak pressures and related temperatures defines a P/T ratio that can be used to characterize the nature of continental subduction responsible for HP-LT metamorphism (Jolivet et al., 2003). Subduction leading to HP-LT metamorphism mostly involves both oceanic and continental material and hence is called here continental subduction (see Jolivet et al., 2003). In the Aegean region, a single northward subduction has been active during the Cenozoic period with different paleogeographic domains (from north to south: Vardar ocean, Pelagonian microcontinent, Pindos ocean, Apulian continent, and the eastern Mediterranean oceanic domain) that were successively subducted and partly accreted to the southern margin of Eurasia (Jolivet and Brun, 2008; Jolivet et al., 2003; Van Hinsbergen et al., 2005). The HP-LT rocks of the Aegean region are taken here as an example to illustrate the linear relationship observed in PT diagrams (Fig. 1). Two metamorphic belts are well documented in the Aegean, from north to south: the Cycladic blueschists formed ca. 40–45 Ma (i.e., grey areas in Fig. 1A), and the Cretan blueschists are Miocene age (25 Ma; black domains in Fig. 1A). In the Cyclades, PT data for Evvia (Ev), Tinos (Ti), Samos (Sam), Syros (Syr), and Syfnos (Syf), respectively, are taken from Shaked et al. (2000), Parra et al. (2002), Will et al. (1998), and Trotet et al. (2001). In the Cretan blueschists, PT data for eastern Crete (EC), central Crete (CC), and western Crete (WC and WCXa, which stands for the Xania site in western Crete) are taken from Theye et al. (1992) and Jolivet et al. (1996). The Cycladic blueschists have yielded scattered ages, ranging from the Late Cretaceous for some eclogites (70 Ma in Syros; Bröcker and Enders, 1999) to the Eocene for most blueschists (Wijbrans and McDougall, 1986, 1988), yielding a mean age for the peak metamorphism at ca. 45 Ma (Jolivet et al., 2003). Younger HP rocks are found farther south in Crete and Peloponese, with a pressure peak dated at 25 Ma in Crete (Jolivet et al., 1996). Readers are referred to Jolivet et al. (2003) for a compilation of PT data and for the ages of peak metamorphism. Peak pressure and related temperature for the Eocene and Miocene blueschists are reported in Figure 1B. PT data from Cycladic blueschists are well aligned within the PT diagram. Two parameters will be used in this paper to characterize this linear relationship: the slope and the temperature at 10 kb of the P/T alignment, respectively, denoted ΔP/ΔT and T@10 kb. The Cycladic alignment yields ΔP/ΔT = 0.056 and T@10 kb = 375 °C. Note that Jolivet et al. (2003) gave a P/T ratio of 0.03 for the Cycladic blueschists because they use lines with an origin at 0 kb and 0 °C, although most of the linear relationships do not pass through the origin of the PT diagram. For the Miocene blueschists, the Cretan alignment yields a higher ΔP/ΔT and a lower T@10 kb.

The causes of the common alignment of peak PT conditions for HP-LT metamorphisms remain poorly understood. Peacock and Wang (1999) used the Japan subduction example to discuss the role of subduction velocity in defining the thermal gradient within the subduction zone, which differs from the classical continental geotherm. Peacock (1993) and van Keken et al. (2002) show that the PT path of the bottom of the subducted oceanic crust is colder for an increasing subduction dip or velocity. Abers et al. (2006) and Roselle and Engi (2002) show more specifically that the thermal gradient within the subduction zone and, thus, the prograde PT path of metamorphic rocks, is a function of the convergence rate, the width of the subduction channel, the location of the rocks within the subduction channel, the rate of shear heating along the subduction zone, and finally, the intensity of convection within the mantle wedge. The PT path for the top of the subducted crust yields very hot conditions due to the efficient heat transfer from the corner flow system on top. In contrast, the PT path for the bottom of the subducted crust yields conditions that are too cold. An intermediate position within the subducted crust yields an accurate estimate of the PT path. Arcay et al. (2007b) also show through 2D numerical simulations that slab surface top temperature is controlled by upper plate thinning processes and hence by mantle wedge flow. However, the quantitative role of these parameters in defining the shape of the prograde PT path and the alignment of peak PT conditions of several HP-LT rocks remains to be quantified. A compilation of PT data from the Mediterranean HP-LT rocks shows a correlation between the P/T ratio and the subduction velocity (Jolivet et al., 2003). Large values for the P/T ratio are likely related to large subduction velocity values. Because this correlation is only based on compilations of natural data, physical explanations of that type of alignment are still lacking and are, therefore, the subject of this paper.

We undertake simple 1D thermal modeling of continental subduction in order to understand the relationship between subduction velocity and the P/T ratio. Previous numerical models commonly compare natural PT paths with modeled PT paths in order to infer the thermal state and burial and exhumation rates from natural data (Bousquet et al., 1997; Burov and Yamato, 2008; Gerya et al., 2002; Huerta et al., 1998; Jamieson and Beaumont, 1988; Jamieson et al., 2004; Roselle and Engi, 2002; Stöckhert and Gerya, 2005; Yamato et al., 2008). The shapes of modeled PT paths are very similar to those of natural rocks but with colder temperatures (Burov and Yamato, 2008; Yamato et al., 2008), highlighting a possible impact of shear heating (Abers et al., 2006; Burg and Gerya, 2005) or radiogenic heat production in the crust (Goffé et al., 2003). Such models successfully demonstrated the necessity of high exhumation rates to explain adiabatic retrograde paths, which are a common feature of HP-LT metamorphic rocks. Moreover, the role of the exhumation velocity is well described by the Peclet number (Pe) for heat transfer (Batt and Braun, 1997; Platt et al., 1998), which is the ratio of the advection time (a function of the exhumation velocity Ve) over the conduction time (a function of the characteristic length). At a high Pe, the retrograde path is mostly adiabatic. Most natural adiabatic retrograde paths of HP-LT metamorphic rocks are, therefore, likely related to high exhumation rates. However, parameters that control the pressure and temperature at peak pressure conditions are little discussed in these previous models. In the present paper, we investigate the role of subduction velocity and dip angle in defining the P/T ratio of HP-LT metamorphic rocks.

One-dimensional (1D) thermal modeling of a thinned continental crust undergoing subduction is first presented where we show the respective roles of subduction velocity and dip angle in defining ΔP/ΔT and T@10 kb values. The model results are then applied to the Aegean example, where the variations through time of the P/T ratio between the Cyclades and the Cretan blueschists are related to variation of subduction velocity and dip angles.

## ONE-DIMENSIONAL MODELING OF HP-LT PROGRADE PT PATH

### Model Setup

We performed 1D modeling of the temperature evolution of a vertical segment within the subducted lithosphere. Here, we only focus on burial history, e.g., the prograde PT path. This 1D model consists of a transient conductive model of a thinned continental lithosphere undergoing subduction. One-dimensional heat conduction is solved in the subducted lithosphere and in a thin diffusive layer on top. This top layer belongs to the overlying plate and accounts for thermal exchange between the subducted lithosphere and the overlying plate, and is explained below. The initial temperature condition is based on a 1D temperature distribution of a passive margin with an age of 120 Ma after breakup, defined by the previous modeling done by Leroy et al. (2008). At that age, the thermal equilibrium is almost reached, and, thus, the lithosphere-asthenosphere boundary is nearly flat. In the present 1D model, we will only study the thinned part of the passive margin, near the continent-ocean boundary. The subducted lithosphere is thus typical of a continental margin with a 10-km-thick continental crust and an overall thickness of 100 km (Fig. 2A). The selected thermal parameters are those used in a previous 2D modeling study (Carry et al., 2009) and are as follows: the capacity is set to 1000 J kg^{−1} K^{−1}, the density in the crust and mantle are, respectively, 2700 kg m^{−3} and 3300 kg m^{−3}, the diffusion coefficients in the crust and mantle are, respectively, 2.1 and 3.0 W m^{−1} K^{−1}, and the radiogenic heat production is set to 10^{−6} W m^{−3} and is only accounted for in the first three kilometers of the thinned continental crust.

The boundary conditions are as follows. A constant mantle heat flux is applied at the bottom, while the top temperature evolves through time. The subduction velocity V_{s} and dip angles, α_{c} and α_{m}, define the time at which the top of the vertical segment within the subducted lithosphere will reach a given burial depth (Fig. 2A). The top temperature at this time is simply the temperature of the overlying lithosphere geotherm at that burial depth. The subduction dip angle is assumed to be 10° (α_{c} = 10°) when subduction occurs in the overlying crust, whereas the dip angle of the subduction in the mantle, α_{m}, is a free parameter and will be systematically changed (Fig. 2A). These boundary conditions are similar to those used by van Keken et al. (2002), with a change in the subduction dip when subducted material enters the overlying mantle. Moreover, the temperature distribution in these previous 2D numerical models is almost unchanged within the overlying lithosphere, where the corner flow is limited at depths lower than 100–150 km. The corner flow, which is mostly controlled by mantle viscosity, has a strong impact in the thermal gradient of the subduction zone and yields a large temperature increase in the top of the subducted crust for pressures higher than 30 kb. In the present study, we do not account for time variation in the overlying geotherm, and, thus, we disregard the role of corner flow. This implies that our model results can only be applied to metamorphic rocks with a peak pressure lower than ~30–40 kb. The steady geotherm of the overlying lithosphere is defined by the basal heat flux q_{m}, the overlying crustal thickness, which is set to a common value for unthinned continental crust at 30 km, and the radiogenic crust thickness, which is one-third of the crustal thickness. The top temperature therefore increases with time and is a function of the subduction velocity and dip angle, V_{s} and α_{m}, respectively. The positions of the top of the 1D subducted lithosphere segment at different times are reported on the 2D pictures of the model setup (Fig. 2A), and the corresponding 1D modeled temperature profiles through the subducted lithosphere are plotted in Figure 2B. This modeling approach is similar to that used by Huerta et al. (1998). However, we do not account for erosion rate in our modeling because we focus here on prograde PT evolution during the early stage of continental subduction, prior to collision and mountain building. At this stage, the convergence between the two plates is mostly accommodated within the subduction zone, where the stacking of HP-LT metamorphic units occurs, hence, leading to low rate or an absence of erosion (Guillot et al., 2003; Hacker et al., 2003a; Hacker et al., 2003b; Jolivet et al., 2005). Starting from a steady state, the temperature profile within the subducted lithosphere evolves through time because of the increase in top temperature, defining classical reversed temperature profiles (Fig. 2B). Note that there is no thermal exchange between the studied structure and the overlying lithosphere because the geotherm of the overlying plate remains constant through time (Fig. 2A). The presence of a diffusive layer on top of the subducted lithosphere allows us to indirectly model the thermal exchange between the overlying lithosphere and the subducted lithosphere. The thickness and the conductivity of this layer, respectively, 3 km and 1.785 W m^{−1} K^{−1}, were adjusted so that the 2D thermal structure of the subduction zone that is interpolated from the 1D profiles at different burial depths is comparable to that of previous 2D numerical models. The advantage of a 1D model is the capacity to quickly run several hundred models in order to systematically analyze the role of subduction velocity and dip angle in defining the P/T ratio of HP-LT rocks. The second advantage of a simple 1D model is having few free parameters, so that it is possible to fully quantify their respective roles. In our problem, we only have three parameters to analyze: (1) subduction velocity V_{s}, (2) subduction dip angle α_{m}, and (3) basal heat flux q_{m}. Finally, note that this 1D model follows the same modeling approach of Carry et al. (2009), in which a detailed justification of all modeling assumptions is given.

### Prograde PT Path

The temperature of a crustal piece initially at 5-km depth within the subducted lithosphere is followed through time, defining the prograde PT path (Fig. 2C). We chose that specific initial depth because most of the metamorphic nappes in modern orogens and, in the Aegean in particular, are made of upper crustal material (Van Hinsbergen et al., 2005). As also discussed earlier, previous numerical modeling has shown that prograde PT paths are strongly controlled by the position of the crustal fragment within the subducted crust (Abers et al., 2006; Roselle and Engi, 2002; van Keken et al., 2002). An intermediate position within the subducted crust was shown to yield an accurate estimate of the PT path, and justified our choice of a crustal piece at 5 km. The particle pressure is simply the lithostatic pressure and, thus, is directly related to the burial depth, which is a function of the subduction velocity and dip angle. The overlying crustal and mantle densities are 2700 kg m^{−3} and 3300 kg m^{−3}, respectively. The prograde PT paths for V_{s} = 1 cm a^{−1}, α_{m} = 30°, and q_{m} = 30 mW m^{−2} are given in Figure 2C. Because the subduction dip angle is lower in the overlying crust than in overlying mantle (α_{c} < α_{m}), the vertical burial velocity in the overlying crust is lower than burial velocity in the overlying mantle. As a consequence, thermal re-equilibration during subduction in the crust is favored, while a more adiabatic evolution for subduction in the mantle is expected. The differences in slopes of the PT path for pressures both lower and greater than 8 kb are thereby explained by the differences in subduction dip angles in the overlying crust and mantle (Fig. 2C). Two parameters can be defined to fully characterize the PT path—the slope ΔP/ΔT, defined for pressure between 10 kb and 20 kb, and the temperature at 10 kb, T@10 kb.

### Role of the Subduction Velocity V_{s}

Figure 3A presents prograde PT paths for the same values of basal heat flux and subduction dip angle (q_{m} = 30 mW m^{−2}; α_{m} = 30°) and for four different values of subduction velocity V_{s}: 0.5 cm a^{−1} (curve 1), 1 cm a^{−1} (curve 2), 1.5 cm a^{−1} (curve 3), and 2 cm a^{−1} (curve 4). The PT path is colder at larger subduction velocities V_{s} because the faster a crustal piece is subducted, the more adiabatic is its transient thermal evolution. Therefore, an increase in subduction velocity decreases the temperature at 10 kb. For example, T@10 kb is 420 °C for V_{s} = 1 cm a^{−1}, while it is only 350 °C for V_{s} = 2 cm a^{−1}. The slope of the PT path ΔP/ΔT is also a function of V_{s}: ΔP/ΔT is 0.04 and 0.05 for V_{s} equal to 1 and 2 cm a^{−1}, respectively. Doubling the subduction velocity thereby leads to a decrease of T@10 kb by 20% and to an increase of ΔP/ΔT by 25%. The subduction velocity V_{s} thus controls both T@10 kb and the slope ΔP/ΔT and, therefore, the overall shape of the P/T path.

### Role of the Subduction Dip Angle α_{m}

Figure 3B presents PT paths for a given value of basal heat flux and for two values of subduction velocity and subduction dip angle: V_{s} = 1 cm a^{−1} and α_{m} = 30° (curve 2, same as curve 2 in Fig. 3A), V_{s} = 1 cm a^{−1} and α_{m} = 50° (curve 2b), V_{s} = 2 cm a^{−1} and α_{m} = 30° (curve 4, same as curve 4 in Fig. 3A), and V_{s} = 2 cm a^{−1} and α_{m} = 50° (curve 4b). At a given velocity (V_{s} = 1 cm a^{−1}), an increase in subduction dip angle from 30° to 50° leads to little decrease in T@10 kb from 420 °C (curve 2) to 400 °C (curve 2b) and to a more significant increase in the slope of the PT path ΔP/ΔT from 0.04 (curve 2) to 0.05 (curve 2b). Because the subduction dip angle is constant when subduction occurs in the overlying crust (α_{c} = 10°), PT paths for the same subduction velocity are similar in the first 8 kb, and, thus, T@10 kb values are very similar, even if the mantle dip angle α_{m} varies (Fig. 3B). At a pressure larger than 8 kb, an increase in subduction dip angle α_{m} leads to an increase in the burial rate and, subsequently, to a more adiabatic evolution of the subducted lithosphere segment. As a consequence, an increase in subduction dip angle leads to an increase in the slope of the PT path ΔP/ΔT.

In summary, the subduction velocity V_{s} controls both T@10 kb and the slope ΔP/ΔT, whereas the subduction dip angle α_{m} only controls the slope ΔP/ΔT. From natural peak PT data, such as those of the Aegean region (Fig. 1B), it would therefore be possible to infer both subduction velocity and subduction dip angle from the P/T ratios. For that purpose, a quantification of the relationship between the subduction kinematics (V_{s} and α_{m}) and the shape of the PT path (T@10 kb and ΔP/ΔT) is required.

## P/T RATIOS AS A FUNCTION OF SUBDUCTION DIP ANGLE AND VELOCITY

### P/T Ratio Diagrams: ΔP/ΔT as a Function of T@10 kb

In order to compile the role of V_{s} and α_{m} in defining the shape of the PT paths, values of ΔP/ΔT have been plotted as a function of T@10 kb in Figure 3C for the set of models in Figures 3A–3B (e.g., α_{m} = 30° and 50° and V_{s}= 0.5, 1, 1.5, and 2 cm a^{−1}). The numbering of the squares in Figure 3C refers to the numbering of PT path curves in Figures 3A–3B. Based on the T@10 kb and ΔP/ΔT values computed by the 1D models, it is possible to draw lines in Figure 3C that are defined by a given subduction velocity (V_{s} isolines are plotted as solid lines) and lines defined by a given subduction dip angle (α_{m} isolines are plotted as dashed lines). In the ΔP/ΔT versus T@10 kb diagram, V_{s} isolines (for V_{s} = 1 and 2 cm a^{−1}) are roughly vertical, whereas α_{m} isolines are oblique. These features highlight the role of V_{s} in defining both T@10 kb and ΔP/ΔT, whereas α_{m} solely controls ΔP/ΔT. To fully describe the role of V_{s} and α_{m} in defining the P/T ratio, we ran 40 models for a given heat flux (q_{m} = 30 mW m^{−2}) in order to construct a ΔP/ΔT versus T@10 kb diagram with five α_{m} isolines (10°, 20°, 30°, 40°, and 50°) and eight V_{s} isolines (0.5, 1, 1.5, 2, 2.5, 3, 3.5, and 4 cm a^{−1}). The corresponding P/T ratio diagram is shown in Figure 4A. The tendency observed in Figure 3C is confirmed with a larger number of numerical results: V_{s} isolines are almost vertical, while α_{m} isolines are oblique. An important feature is that V_{s} isolines are distinct, even with differences in V_{s} as low as 0.5 cm a^{−1}. More specifically, for α_{m} = 30° (30° dashed isoline), T@10 kb progressively decreases with an increasing V_{s}: 475, 415, 380, 350, 330, 315, 300, and 285 °C, respectively, for V_{s} = 0.5, 1, 1.5, 2, 2.5, 3, 3.5, and 4 cm a^{−1}. Differences in T@10 kb are more important for low V_{s} values than for large V_{s} values. This feature is explained as follows. At low V_{s} (0.5 cm a^{−1}), thermal evolution of the subducted lithosphere is close to a steady state. Therefore, a small increase in V_{s} will strongly change the transient evolution of the subducted lithosphere, explaining a significant decrease of T@10 kb for a small increase in V_{s}. For a large V_{s} (>2 cm a^{−1}) value, the burial rate is high, and, thus, the burial time is low compared to the characteristic heat conduction time, explaining a quasi-adiabatic evolution of the subducted lithosphere. Thus, an increase in V_{s} will not change this transient thermal evolution because the asymptotic adiabatic evolution is almost reached.

### Role of Basal Heat Flux

The role of q_{m} in defining the shape of the PT path can now be discussed through the P/T ratio diagram. Figures 4B and 4C, respectively, show P/T ratio diagrams for lower (25 mW m^{−2}) and higher (35 mW m^{−2}) basal heat flux compared to 30 mW m^{−2} (the reference value discussed above; Fig. 4A). A decrease in basal heat flux leads to (1) a shift toward a lower temperature for V_{s} isolines and (2) a shift toward higher ΔP/ΔT values for α_{m} isolines (Fig. 4B). A decrease in the basal heat flux q_{m} indeed leads to a decrease in the overlying plate geotherm and to a decrease in the initial geotherm of the subducted lithosphere. As a consequence, a decrease in q_{m} leads to a decrease in the temperature gradient within the subducted lithosphere and, thus, a shift toward a lower T@10 kb for the V_{s} isolines and a shift toward larger values of ΔPΔT for the α_{m} isolines (Fig. 4B). The opposite is observed for an increase in q_{m}. An important feature is that the changes of the position for both V_{s} and α_{m} isolines are large, even for a decrease and/or increase of a basal heat flux of 5 mW m^{−2}.

## APPLICATIONTOTHE AEGEAN BLUESCHISTS

### Method to Infer Subduction Velocity and Dip from Natural PT Data

Figure 5 schematically explains the method used in this paper to infer subduction kinematics at the time of the peak metamorphism from natural thermobarometric data. Natural continental subduction implies HP-LT metamorphism of several units. These units are assumed to belong to the same slab, if they have the same ages for peak pressure conditions. If natural thermobarometric data are aligned in the PT diagram, they yield a P/T ratio (i.e., ΔP/ΔT and T@10 kb values) for the considered HP-LT units. For a given value of the basal heat flux q_{m}, the plot of these ΔP/ΔT and T@10 kb values in the P/T ratio diagram permits one to infer both the subduction velocity and subduction dip angle. If the value of the subduction velocity V_{s} is consistent with an independent estimate of convergence velocity, from a paleogeographic reconstruction, for example, it validates the assumed value of the basal heat flux q_{m}. In other cases, either geodynamic reasons (e.g., subduction rollback, oblique subduction, etc.) or variation in basal heat flux through time can explain the differences between the subduction velocity V_{s} and the convergence velocity V_{c}. A decrease in q_{m} leads to a shift of V_{s} isolines to a lower T@10 kb (Fig. 4), and, thus, the same natural values of ΔP/ΔT and T@10 kb can yield lower values of V_{s} for a lower q_{m}. The Aegean example is now used to illustrate this procedure (Fig. 6).

### Application to the Eocene and Miocene Aegean Blueschists

The Aegean blueschists were first classified with respect to the ages of peak metamorphism—Eocene Cycladic blueschists and Miocene Cretan blueschists. As discussed in the introduction, these two blueschist belts yield well-defined trends within the PT diagram (Fig. 1) with different ΔP/ΔT and T@10 kb values. Cycladic blueschists yield a ΔP/ΔT of 0.056 (±0.02) and a T@10 kb of 375 °C (±50 °C), while the Cretan blueschists give a much larger ΔP/ΔT of 0.08 (±0.015) with a lower T@10 kb at 315 °C (±20 °C, Fig. 1). Errors in ΔP/ΔT and T@10 kb are a result of errors in the PT estimates. These ΔP/ΔT and T@10 kb values for both Cycladic and Cretan blueschists have been plotted on the P/T ratio for q_{m} = 30 mW m^{−2} in Figure 6A. The Cycladic trend yields V_{s} of 1.5 cm a^{−1} (±1 cm a^{−1}) and a dip angle of 50° (±30°). The large errors in V_{s} and α_{m} estimates are due to errors in the PT estimates. This estimate of V_{s} is very similar to the estimates of the horizontal Africa-Europe convergence velocity, which is ~1–1.5 cm a^{−1} (Dewey et al., 1989; Rosenbaum et al., 2002). Therefore, this good consistency indicates that (1) the simple 1D model used here is relevant to infer the kinematics of continental subduction at the time of peak metamorphism from natural PT data, and (2) the selected value of the basal heat flux is appropriate to describe the thermal regime of the Aegean subduction. Moreover, the use of different values of q_{m} (25 mW m^{−2} or 35 mW m^{−2}, Fig. 4) yields two very different V_{s} values (0.7 cm a^{−1} or 2.5 cm a^{−1}, respectively), which are inconsistent with the known values of the Africa-Europe convergence rate. The same procedure is applied for the Cretan blueschists with the same basal heat flux (q_{m} = 30 mW m^{−2}), yielding a V_{s} of 2.75 cm a^{−1} (±0.5 cm a^{−1}) and a subduction dip angle greater than 50° (Fig. 6A). This V_{s} value is very different from the value inferred for the Cyclades, even if we account for the large errors in the V_{s} estimates.

The model results therefore imply that the differences in the P/T ratio between Crete and Cyclades are primarily due to large differences in subduction velocity, which was greater during the Miocene subduction than during the Eocene subduction. These findings are in good agreement with the geodynamic evolution of the Aegean from the Eocene to Miocene, as discussed below.

### Implication for the Aegean Slab Dynamics through Time

Figure 6B summarizes the model results in terms of subduction velocity and dip angle. The Eocene Cycladic subduction occurred at a rate V_{s} of 1.5 cm a^{−1}, which is comparable to the Africa-Europe convergence rate V_{c} of 1–1.5 cm a^{−1} (Dewey et al., 1989; Rosenbaum et al., 2002). This indicates a continental subduction with low or no rollback, where the subduction velocity is almost equal to the convergence velocity. The model results indicate contrasting features for the Miocene Cretan subduction that occurred at a much larger rate (V_{s} ~2.75 cm a ^{−1}). This larger subduction velocity can reflect (1) a larger convergence rate V_{c} and/or (2) a component of rollback velocity V_{rb} that is added to the convergence rate V_{c} to define the subduction velocity V_{s} = V_{rb} + V_{c}. Previous studies that discussed the rate of convergence between Africa and Europe do not show an increase of that rate between the Eocene and Miocene period (Dewey et al., 1989; Rosenbaum et al., 2002). The last hypothesis thus seems to be the most appropriate to explain the increase in subduction rate for the Miocene subduction—a component of rollback that has to be of the order of V_{rb} = 1.25 cm a^{−1}. Previous studies have estimated the amount of rollback to be ~300 km from the Cyclades to the Cretan blueschists between 45 Ma and 25 Ma (e.g., a duration of 20 Ma), yielding 1.5 cm a^{−1} of rollback velocity (Jolivet and Faccenna, 2000). Recent estimates made by Jolivet and Brun (2008) yield 200 km of Aegean slab rollback between 75 and 35 Ma, and 500 km between 35 and 0 Ma. This implies a rollback rate of 0.5 cm a^{−1} during the subduction of Cycladic blueschist and 1.5 cm a^{−1} during the subduction of Cretan blueschist. Using 1 cm a^{−1} for the Africa-Europe convergence rate V_{c} (Dewey et al., 1989; Rosenbaum et al., 2002), these estimates yield consistent values of subduction velocities V_{s} at ~1.5 cm a^{−1} for the Cylcadic subduction and 2.5 cm a^{−1} for the Cretan subduction.

The consistency between the subduction velocity and dip angle, inferred from the trend defined by natural PT data and the independently known variation through time of the Aegean slab kinematics, validates the use of simple 1D thermal models of prograde PT evolution of HP-LT rocks to infer kinematics of continental subduction.

## CONCLUSIONS

(1) One-dimensional thermal modeling of continental subduction allows us to infer the subduction velocity and dip angle at the time of peak-pressure metamorphism from natural PT data of HP-LT metamorphic units. This new method applies as follows. First, a compilation of thermobarometric data for peak-pressure metamorphism of HP-LT units of the same age must be plotted in a PT diagram in order to infer a possible trend underlined by these PT data. This trend would be defined by a slope, denoted ΔP/ΔT, and an origin defined at 10 kb, denoted T@10 kb. Second, the 1D thermal modeling we performed permits the quantification of the role of both subduction velocity and dip angle in defining the slope and the origin of a PT path, for a given value of the basal heat flux. Third, the comparison of natural data (ΔP/ΔT and T@10 kb) and the model results permits us to directly infer both the subduction velocity and dip angle that prevailed at the time of the peak-pressure metamorphism.

(2) This new approach is validated in the Aegean region for two metamorphic belts of Eocene and Miocene age, respectively, the Cyclades and Crete. The differences in P/T ratios between Crete and Cyclades underlined by the HP-LT rocks have been primarily interpreted in terms of variation in subduction velocity through time. The differences in subduction velocity are likely related to a rollback velocity of the order of 1.25 cm a^{−1}. These findings are consistent with previous estimates obtained by independent methodology.

(3) Natural peak-pressure thermobarometric data can thus be used as a proxy to quantify subduction velocity and dip angle at the time of the peak metamorphism. This new method can be useful to discuss spatial and temporal changes of the continental subduction and, therefore, can give further constraints on 3D geodynamic evolution through time of modern orogens.

These 1D models provide a first step toward a broader modeling of prograde and retrograde PT of HP-LT metamorphic units. The application of this new approach would be useful in a fully coupled, 2D, thermomechanical modeling of continental subduction (Arcay et al., 2007a; Burov and Yamato, 2008; Gerya et al., 2002; van Keken et al., 2002; Yamato et al., 2008), in which one can account for the shear heating, frictional heating, and heat transfer between the subducted lithosphere and the overlying lithosphere.

The authors would like to thank Laurent Jolivet for his constructive comments on an early version of the manuscript. Jean Pierre Brun and Pierre Gautier are thanked for fruitful discussions on the geodynamic evolution of the Aegean region. Detailed review of Romain Bousquet and comments of Raymond M. Russo, *Lithosphere* science editor, helped us to revise the manuscript.