Seismic visibility of a deep subduction channel – insights from numerical simulation of high-frequency seismic waves emitted from intermediate depth earthquakes

Return flow in a deep subduction channel (DSC) has been proposed to explain rapid exhumation of high pressure–low temperature metamorphic rocks, entirely based on the fossil rock record. Supported by thermo-mechanical models, the DSC is envisioned as a thin layer on top of the subducted plate reaching down to minimum depths of about 150 km. We perform numerical simulations of highfrequency seismic wave propagation (1–5 Hz) to explore potential seismological evidence for the in situ existence of a DSC. Motivated by field observations, for modeling purposes we assume a simple block-in-matrix (BIM) structure with eclogitic blocks floating in a serpentinite matrix. Homogenization calculations for BIM structures demonstrate that effective seismic velocities in such composites are lower than in the surrounding oceanic crust and mantle, with nearly constant values along the entire length of the DSC. Synthetic seismograms for receivers at the surface computed for intermediate depth earthquakes in the subducted oceanic crust for models with and without DSC turn out to be markedly influenced by its presence or absence. While for both modelsP andS waveforms are dominated by delayed high-amplitude guided waves, models with DSC exhibit a very different pattern of seismic arrivals compared to models without DSC. The main reason for the difference is the greater length and width of the low-velocity channel when a DSC is present. Seismic velocity heterogeneity within the DSC or oceanic crust is of minor importance. The characteristic patterns allow for definition of typical signatures by which models with and without DSC may be discriminated. The signatures stably recur in slightly modified form for earthquakes at different depths inside subducted oceanic crust. Available seismological data from intermediate depth earthquakes recorded in the forearc of the Hellenic subduction zone exhibit similar multi-arrival waveforms as observed in the synthetic seismograms for models with DSC. According to our results, observation of intermediate depth earthquakes along a profile across the forearc may allow to test the hypothesis of a DSC and to identify situations where such processes could be active today.


Introduction
Since the establishment of plate tectonics, most ideas on structure and properties of the plate interface along subduction zones were inspired by textbook cartoons, showing a giant thrust fault with its position being marked by Benioff zone seismicity.The narrow fault or shear zone concept also underlies thermal models of subduction zones (e.g.Turcotte and Schubert, 2002;Van Den Beukel and Wortel, 1988;Peacock, 1996;Syracuse et al., 2010), with a downward single pass flow beneath a rigid mantle wedge in the forearc.This starting model became challenged by the finding that high pressure over temperature (P /T ) metamorphic rocks are very widespread within continental crust (e.g.Ernst and Liou, 1999;Tsujimori et al., 2006;Agard et al., 2009;Guillot et al., 2009;Ota and Kaneko, 2010), obviously having been exhumed from a subduction zone.While in many cases some cooling associated with decompression requires Fig. 1: Envisaged geometry of deep subduction channel (DSC), including subduction erosion and underplating (after Stöckhert, 2002).Insets show scale-invariant block in matrix (BIM) structure and typical pressure-temperature (P-T) path derived for rapidly exhumed high P /T metamorphic slices.exhumation as relatively small bodies with heat transfer into a cooler environment (e.g.Schertl et al., 1991;Stöckhert et al., 1997;Reinecke, 1998;Groppo et al., 2007;Groppo and Castelli, 2010;Krebs et al., 2011;Angiboust et al., 2012), geochronological results pose narrow time windows for typical timescales of exhumation, asking for displacement rates comparable to plate velocity (e.g.Gebauer et al., 1997;Rubatto and Hermann, 2001;Baldwin et al., 2004;Little et al., 2011).Rapid exhumation with concomitant cooling reconciles neither with uplift and erosion nor with crustal extension (Platt, 1993;Ring and Brandon, 1999).Moreover, geochronological results suggest that in some cases a large part of exhumation has taken place during ongoing subduction, prior to collision (e.g.Lardeaux et al., 2001;Rubatto et al., 2011).These observations have led to postulation of a deep subduction channel (DSC) with forced return flow as a possible mechanism for exhumation (Fig. 1).
The concept of a subduction channel was first applied to explain structure and P -T paths of the high P /T metamorphic Franciscan mélange in California by Cloos (1982), and elaborated by Shreve and Cloos (1986) and Cloos and Shreve (1988a, b).It was inferred to play a role in the frontal forearc beneath an accretionary wedge, competing with models based on critical wedge geometry (e.g.Platt, 1993).Identification of ultra-high-pressure (UHP) metamorphic rocks starting in the 1980s (e.g.Chopin, 1984;Schertl et al., 1991;Schreyer, 1995;Harley and Carswell, 1995), and in particular the finding that these rocks undergo very rapid exhumation (e.g.Gebauer et al., 1997;Rubatto and Hermann, 2001;Baldwin et al., 2004;Parrish et al., 2006), inspired tectonic models invoking a DSC reaching to greater depth.Forced return flow of material from depths exceeding 100 km, as schematically shown in Fig. 1, seemed to provide a feasible alternative for generation and exhumation of smaller UHP metamorphic slices, distinct from coherent crustal or lithospheric sec-tions of the downgoing plate during collision (e.g.Chemenda et al., 1995;Hacker et al., 2000;Kylander-Clark et al., 2012).Numerical simulations (e.g.Burov et al., 2001;Gerya et al., 2002;Warren et al., 2008) generating bulk flow patterns and P -T paths for individual material points, grossly consistent with natural record and structures of orogenic belts (e.g.Engi et al., 2001;Stöckhert and Gerya, 2005), support the principal feasibility of the DSC concept.
Geometry and width of a DSC remain poorly constrained by natural observations (e.g.Warren, 2013).A lower bound to thickness is given by the size of high P /T metamorphic slices supposed to be carried by forced return flow in a lowviscosity matrix, and an upper bound by the requirement of high rate of return flow (e.g.Gerya and Stöckhert, 2006).These bounds are assumed to be on the order of 1 and 10 km, respectively.The fact that high P /T metamorphic rocks can undergo cooling during extremely rapid exhumation requires slices or blocks of limited size, allowing effective heat transfer into a cooler environment.Embedding into a cooler environment, by whatever mechanism and kinematics, requires large amounts and rates of displacement, which in turn imply high strain rates and a complex structural evolution.A DSC is envisaged to obey these conditions, invoking chaotic flow in a heterogeneous rock assemblage, and effective mixing of material of different provenance on different length scales (e.g.Gerya andStöckhert, 2002, 2006).High P /T metamorphic tectonic mélanges are therefore likely candidates to represent material properties in a fossil subduction channel.
When viewed on mesoscopic scale, mélanges frequently exhibit block-in-matrix (BIM) structures.At best slightly deformed blocks of variable provenance are embedded in a highly strained matrix, which has undergone prolonged deformation (Bebout, 2007;Festa et al., 2010;Grigull et al., 2012).The evolution of this spectacular structure is still poorly understood.Importance of sedimentary processes at the surface versus tectonic processes at depth, or both in combination, remains a matter of debate (Dilek et al., 2012).In any case, BIM structures imply contrasts in mechanical properties -the blocks undergoing no notable permanent deformation while the matrix must deform at very high strain rate, despite identical conditions in terms of temperature and stress.The mechanism allowing for deformation of finegrained polyphase rocks at very low differential stress and high strain rates is dissolution precipitation creep (Wassmann and Stöckhert, 2012;Grigull et al., 2012).Consequently, mélanges with BIM structures reflect low viscosity, which may be an important factor to enable forced return flow in a DSC.
Typical bulk material properties in a subduction channel at an arbitrary depth can hardly be predicted.Exhumed high P /T metamorphic mélanges, inferred to represent a fossil DSC, demonstrate that a variety of rock types, both in terms of bulk composition and metamorphic mineral assemblage, occur intermingled on various length scales (Krebs et al., 2011;Grigull et al., 2012), ranging from map scale (e.g Angiboust et al., 2012) down to the scale of the individual outcrop (e.g.Federico et al., 2007).As basic eclogite and serpentinite probably represent extremes, the amplitude of heterogeneity is constrained by the contrast in mechanical properties between these rock types.Therefore, we assume a composite material with eclogite blocks embedded in a serpentinite matrix to account for heterogeneity in the DSC, with a block size roughly 1 order of magnitude smaller than the wavelength of P waves.Taking these extremes, other heterogeneities are assumed to produce a much weaker effect in terms of seismic velocity.
Even for a uniform and specified bulk composition, it would not be possible to predict physical rock properties, as adaption of phase assemblages to changing temperatures and pressures is controlled by kinetics and generally delayed both during burial and exhumation.Blueschists, eclogites, or in particular ultrahigh pressure metamorphic rocks at the surface demonstrate substantial disequilibrium.As such, material in a DSC must be highly heterogeneous, both in bulk composition and in mineral phase assemblage.
The DSC concept was developed to explain the rock record (i.e.structural aspects of high P /T metamorphic terranes, petrologic observations, and geochronological constraints).Whether forced return flow is commonplace in subduction zones, or whether it is a special case restricted to certain stages or geometrical properties, so far remains an unresolved question.While structural geology and metamorphic petrology deal with the product of long-term processes taking millions of years, seismological observations may provide evidence for the presence of a currently active DSC in situ.Here, we investigate by numerical simulations whether a DSC can be detected using high-frequency seismic waves (up to 5 Hz) recorded at the surface.
Guided waves in a subduction zone have previously been simulated by Martin et al. (2003), Martin et al. (2005) and Martin and Rietbrock (2006) using finite-difference modeling.They studied the influence of slab geometry and earthquake locations on the properties of guided waves but their model did not consider the presence of a DSC.Our simulations extend previous work by Essen et al. (2009) who investigated the influence of a DSC on comparatively lowfrequency (< 1 Hz) guided waves and assumed constant seismic velocities in predefined subregions of the model.As effect of a DSC, the authors identified a substantial shift of the maximum amplitude of guided waves from the area where the oceanic crust enters the mantle towards an area in the forearc located above the mantle wedge.Owing to the lowfrequency signals, differences in the internal structure of the guided wave signals could not be analyzed.
The geometry of our simple seismic model is based on that of the Hellenic subduction zone, for which geological and seismological phenomena (Meier et al., 2007) are tentatively ascribed to anomalous material extruded from a DSC beneath Crete.Seismic velocities are calculated based on pre-dictions on temperature and pressure field and assumptions on rock composition and mineral assemblage.For the subduction channel, a BIM structure is assumed composed of a serpentinite matrix and eclogite blocks, for simplicity.In the simulations, we consider seismic waves recorded in the forearc and emitted from intermediate depth earthquakes in the subducted oceanic crust or the lithospheric mantle below.These waves, spending the major part of their travel path in subducted lithosphere and overlying mantle wedge, are expected to be affected by presence of an intervening DSC with distinct properties.For comparison, we also show examples of observed waveforms from intermediate depth earthquakes recorded in the forearc of the Hellenic subduction zone.

A seismic velocity model of a subduction zone
Seismic velocities of mantle rocks are controlled by the physical properties of their mineralogical constituents which in turn depend on temperature and stress state.Due to phase transitions during subduction which may be delayed by kinetic constraints on the mineral reactions, physical rock properties may also depend on the dynamics of the subduction process.
While seismic velocities in subducted oceanic crust or within an inferred DSC are subject to great uncertainty, typical seismic velocities in both subducted mantle lithosphere and mantle wedge are fairly well known from seismic tomography.In order to obtain consistent, pressure and temperature dependent values of P velocities, S velocities and density in our model, we chose to calculate these values for presumed standard rock compositions.

Temperature model
As a first step towards seismic velocities, a temperature model of the Hellenic subduction zone is set up.We consider a vertical cross section running in north-south direction from the Cyclades to south of Crete.The structure along this transect is quite well known from previous seismological studies (Meier et al., 2007;Endrun et al., 2004;Snopek et al., 2007).Basic building blocks are from south to north (Fig 2,top): a 90 km thick continental African lithosphere in the south with 32 km thick crust on top, an accretionary wedge underlain by subducting oceanic crust and oceanic mantle lithosphere dipping at an angle of 15 degrees down to a depth of 60 km and then bending to an angle of 30 degrees at greater depths, and 60 km thick Aegean lithosphere with 20 km thick crust in the north.The subducted oceanic crust is assumed to be 7 km thick.
Based on this structural model a temperature model was developed using analytic expressions for temperature in different regions of the subduction zone.Although it is not as accurate as a numerical fluid dynamical model it approaches such models very well as shown by Davies (1999).To calculate the complete temperature field using analytical expressions, the subduction zone is divided into several regions (Fig 2,top) -regular crust and lithosphere of overriding and subducted plate (foreland and overriding lithosphere): use of a parabolic temperature-depth relation based on static heat conduction with internal heating down to a depth of 24 km and a temperature at the base of the lithosphere of 1250 • C; -the mantle below the subducting plate as well as the mantle wedge below the overriding Aegean lithosphere and above the slab: temperature increase following an adiabatic gradient of 0.4 K km −1 ; -the forearc zone including the accretionary wedge and the upper part of the slab dipping at an angle of 15 • (slab 1): here we use analytic expressions derived by Royden (1993) for the calculation of temperature in eroding orogenic belts and accretionary prisms with prescribed erosion and accretion rates and a velocity of the ingoing lithosphere of 40 mm yr −1 ; -the region inside the more steeply dipping slab including a boundary layer surrounding the slab to allow for diffusion of heat (slab 2): here we assume timedependent heat conduction perpendicular to the slab dip and advection of temperature parallel to slab dip.An analytical solution to this problem has been given by Davies (1999); -a region between the forearc zone and the steeper part of the slab where none of the above expressions are valid (gap): here we interpolate temperature using neighboring values.
With typical values for thermal diffusivity of 10 −6 m 2 s −1 and thermal conductivity of 2.5 W/(K m), we obtain the temperature field shown in Fig. 2, bottom.

Predicting seismic velocities from mineral properties
We follow an approach by Stixrude and Lithgow-Bertelloni (2005) (SLB2005).They use an expression of the mineral's free energy consisting of two parts: a "cold" part which contains an expansion of free energy to fourth order in the finite strain tensor and does not depend on temperature, and a "thermal" part which specifies the free energy of the vibrations of the crystal lattice in the quasi-harmonic approximation.Besides its explicit temperature dependency the thermal part also varies with strain via the lattice vibrational frequencies.The formulation is completely anisotropic and allows for the calculation of the fully anisotropic tensor of elastic constants if corresponding mineral data are available.The first derivative of free energy with respect to a small change of finite strain due to the passage of seismic waves yields the stress tensor and a further derivative gives the tensor of elastic constants.By evaluating the expressions for the elastic constants at conditions where measurements are available (zero finite strain), the expansion coefficients of the free energy can be determined.For the current application, we assume isotropic finite strains and calculate the equivalent isotropic elastic constants of the medium expressed as the bulk (K) and shear modulus (G).To quantify the thermal part of free energy, SLB2005 assume a parabolic Debye density of vibrational states with a maximum frequency that determines the Debye temperature D .This permits an explicit representation of the internal energy, entropy and specific heat of a mineral by integrals which can be numerically evaluated.Strain dependency of the thermal part comes in via the Grüneisen tensor which is the derivative of the lattice vibrational frequencies with respect to strain.For isotropic strains we get the well-known Grüneisen constant.A complete description also requires strain derivatives of the Grüneisen tensor which possesses the same symmetries as the elasticity tensor.In the isotropic case it can be represented by two constants (similar to K and G) which we denote by q and η S .In summary, the SLB2005 theory requires the specification of the following nine parameters at ambient conditions: molar volume V , bulk modulus K and its pressure derivative dK / dP , shear modulus G and pressure derivative dG / dP , Debye temperature D , Grüneisen constant γ , Grüneisen strain derivatives q and η S .
Calculating anisotropic elastic properties from mineral properties as described above using SLB2005 theory is possible but requires much more detailed data on mineral properties.For the cold part, the full tensor of elastic constants and their pressure derivatives at ambient conditions are needed.For the thermal part, data for the Grüneisen tensor and its strain derivatives at ambient conditions are required.Since anisotropy is typically small, a more practical procedure would be to derive isotropic elastic constants from mineral properties and then add anisotropic perturbations with the desired symmetry properties.
A mineral physics approach to calculate seismic velocities for a subduction zone model may be regarded as of limited value, given the uncertainties in terms of mineral assemblage.It should be noted, however, that these uncertainties mainly affect the subducted oceanic crust, a phase diagram for MORB being depicted in Fig. 3.At any rate, effects related to sluggish kinetics in metamorphic reactions or temperature gradients are expected to be much more prominent compared to those related to uncertainties in the mineral database.Low velocity of oceanic crust compared to surrounding mantle as well as the steady increase of velocities with metamorphic grade leading to reduction of the velocity contrast with depth, are essential for the present modeling of seismic wave propagation.

Mineral database
Mineral physics data are taken from the very comprehensive database collected by Holland and Powell (1998) (HP98), which, however, does not directly list the nine parameters needed by the SLB2005 theory.Data associated with shear deformations are completely missing and instead of Debye temperature, Grüneisen constant and q, it provides the entropy at ambient conditions, the thermal expansivity α and the specific heat c P as a function of temperature at atmospheric pressure.Since there is a strong dependency of α(T ) and c P (T ) on the Grüneisen constant γ and its bulk strain derivative q, we applied a grid search to find those values of γ and q that best reproduced thermal expansivity and specific heat given in the HP98-database.Similarly, Debye temperature was determined by a fit to the given entropy value in HP98.This procedure worked well for most minerals (with few exceptions) with errors of less than 5 percent.Finally, to obtain values for the parameter η S associated with the shear part of the Grüneisen strain derivative tensor, we again performed a grid search for that value of η S which best reproduces the quantity 1/(−Gα)(∂G/∂T ) given in tables by Hacker et al. (2003b).In this way, a database of mineral parameters was compiled from which physical properties of subduction zone minerals were calculated using the SLB2005 theory.

Phase diagrams
A major difficulty of estimating physical properties of mantle rocks is the uncertainty about which mineral assemblages (phases) are stable at given temperature and pressure.Usually, it is assumed that only thermodynamically stable phases are present.In that case, it is sufficient to determine phase diagrams for the relevant rock facies.However, even this step is plagued by uncertainties because the rock compositions are not known very well.Even worse, observations of ultra-highpressure rocks at the surface clearly testify that metastable phases may persist in earth's mantle for a long time and that the hypothesis of thermodynamic stability is violated.This is particularly true for subduction zones where exceptionally low temperatures inside the slab may strongly hamper phase transformations as the material moves to greater depths.For simplicity, we chose to determine the stable phase using phase diagrams calculated on the assumption of thermodynamic equilibrium.Hacker et al. (2003b) presents phase diagrams including volume fractions of the minerals contained in the phase for a typical mid-ocean ridge basalt composition, which we apply to rocks in the subducted oceanic crust, and for a harzburgitic composition which we use for mantle rocks (Fig. 3).We do not perform calculations for rocks in the shallow crust since usually density and seismic velocities are fairly well known there from seismic experiments.
It should be noted that the phase diagrams by Hacker et al. (2003b) assume MORB and harzburgite to be 100 percent water saturated.This leads to the appearance of many hydrated phases in the lower left part of the phase diagrams.To avoid these phases appearing in the subducting mantle part, we instead use dry harzburgite in this region of the model.Dry harzburgite only exhibits the phases F, G and K of the Harzburgite phase diagram shown in Fig. 3.
The geological record of exhumed high P /T metamorphic rocks (Stöckhert, 2002) indicates that a DSC could reach down to a depth of at least 130 km.However, since the existence of a DSC is linked to the release of water from the underlying oceanic crust, a consistent petrological model should allow the presence of some water-bearing minerals at these depths.The phase diagram by Hacker et al. (2003b) predicts the transition from the water-bearing zoisite eclogite to dry diamond eclogite at a pressure of about 3.3 GPa for temperatures between 500 and 700 • C.This pressure corresponds to a depth of about 105 km which is much less than required for a deep-reaching subduction channel.Assuming a kinetic delay of the transition owing to the low temperatures in the oceanic crust, we shifted the zoisite-diamond eclogite transition in the MORB phase diagram to a pressure of 4.7 GPa or equivalent depth of 150 km.

Evaluation of seismic wave velocities and density
Seismic velocities are determined on a dense grid that covers the cross section of the Hellenic subduction zone shown in  Temperature at the grid nodes is calculated as described above (see Fig. 2).Depending on the bulk compositions at a given grid point (continental crust, MORB, harzburgite, dry harzburgite), we look up the stable phase and its modal composition in the corresponding phase diagram (Fig. 3) and calculate elastic moduli and density for each constituting mineral using the theory of SLB2005 and the corresponding parameters in the mineral database.Effective elastic moduli of the phase are obtained by a volume fraction weighted average of the individual reciprocal elastic moduli (Reuss average).Density is calculated by averaging the individual density values according to the volume fractions.This procedure is certainly not satisfactory as it is known that the effective elastic moduli may have varying values depending on the spatial distribution of the modal constituents inside the phase.We have chosen to display here the Reuss averaged values to give an impression of the variations of seismic wave speeds inside a subduction zone.Hashin-Shtrikman bounds (Hashin and Shtrikman, 1963) can be readily calculated and provide an estimate of the variability of these values.Since, however, seismic wave speeds derived from tomographic imaging represent averages over rock volumes of kilometer size, a homogenization scheme is needed to obtain accurate estimates for the effective elastic moduli of such a rock volume.This is particularly true in subduction zones where large changes of temperature and also stress may occur within small distances.Clues of the spatial distribution of rocks inside a seismically "visible" rock volume may be taken from numerical models or geological outcrops.
A weak but ubiquitous anisotropy of the oceanic mantle lithosphere is observed in seismological studies represented by a 3 % radial anisotropy in standard earth model PREM (Dziewonski and Anderson, 1981), implying anisotropy of the mantle part of the subducted slab.Also, deformation in the mantle wedge may create anisotropy.In this model, anisotropy is neglected to avoid complications in the seismic wave field that could be erroneously attributed to the presence of a DSC.Moreover, anisotropy of slab and mantle wedge will not notably affect guided waves, in contrast to anisotropy in a DSC or subducted oceanic crust.
In Fig. 4 we show the resulting distribution of compressional and shear wave speed as well as density and quality factor down to a depth of 150 km.Velocities vary from values below 6.0/3.0 km s −1 mostly in the crust to values of 8.7/4.8 km s −1 in the fastest parts of the slab, respectively.We can clearly recognize the gradual velocity increase in the oceanic crust due to release of water.In the wedge above the oceanic crust we obtain very low wave speeds caused by the presence of hydrated phases.They are partially slower than in the oceanic crust leading to a positive jump of velocity at the plate contact down to a depth of 40 km.For density, we obtain a picture that correlates very well with the seismic wave speeds.Hydrous phases now show up as regions of lower density.At about 80 km depth the oceanic crust becomes increasingly dense due to the appearance of various eclogite phases.For some numerical calculations of wave propagation shown later, we also accounted for attenuation.Attenuation in the mantle primarily depends on temperature, pressure and water content and the activation energies and activation volumes of the dissipating mechanisms.We use the following formula to calculate the inverse quality factor for our model (Karato and Jung, 1998): where E D and V D are the activation energy and volume for dissipation of dry rock, respectively, and E W and V W are their counterparts for water-induced dissipation mechanisms.
C OH is the water content in ppm H/Si, r the water content exponent and A D and A W , respectively, are factors determining the contributions of the two terms to attenuation.α is the frequency exponent and A 0 is a pre-exponential factor chosen in a way to produce seismologically reasonable values of the quality factor.Values chosen for the various constants are listed in Table 1.The quality factor is generally high in the area of the subducted slab due to the low temperatures prevailing there (Fig. 4).For this reason, we do not expect strong effects of attenuation on the wave propagation characteristics in the slab.

Effective seismic velocities in a deep subduction channel
Seismic velocities in the DSC require separate considerations.As mentioned in the introduction, we make the simplifying but reasonable assumption that the material inside the DSC is composed of eclogite blocks of various sizes floating in a low-viscosity serpentinite matrix (block in matrix structure).Since seismic waves recorded at the surface have wavelengths of a few kilometers and the typical size of the blocks is about one order of magnitude smaller, it is suffi-cient to consider effective seismic velocities of the composite material.They are computed using a theory of Sabina and Willis (1988) which is applicable up to frequencies where the effective wavelength approaches the size of the blocks.The theory delivers complex values for the isotropic elastic constants of an equivalent effective medium from which frequency-dependent values for wave velocities and scattering attenuation can be calculated.Controlling input parameters for the theory are the elastic constants and density of the constituting materials eclogite and serpentinite and the frequency-size distribution of the blocks.For eclogite, we used a P wave velocity of 8.36 km s −1 , an S wave velocity of 4.67 km s −1 and a density of 3640 kg m −3 .P and S velocities as well as density of serpentinite were taken from  2010) is significantly higher.Moreover, antigorite and deformed serpentinite exhibit strong anisotropy (Bezacier et al., 2010) which is not accounted for in our simulations.Since samples taken from outcrops of exhumed fossil DSCs indicate a complex flow of the serpentinite matrix around eclogite blocks, we argue that the large-scale anisotropy of the BIM structure is much smaller than the one measured on pure samples in the laboratory.Nevertheless, remaining anisotropy may be a further important indicator of the existence of a DSC.
We have done several calculations using different frequency-size distributions of the blocks that were determined by Grigull et al. (2012) from outcrops of exhumed fossil DSCs.The results show that the seismic velocities of a BIM compound lie between the end-member values for blocks and matrix and mainly depend on the relative amount of blocks and their size distribution.Since the velocity of serpentinite is relatively low and the fraction of eclogite blocks small, the effective velocity of the DSC is found to be less than the velocities of the surrounding mantle and oceanic crust.In addition, petrophysical calculations suggest that seismic velocities of serpentinite matrix are nearly constant over the entire length of the DSC because the effects of temperature and pressure approximately cancel each other.All these considerations lead to a DSC model with P wave speed of around 6.6 km s −1 and S wave speed of 3.3 km s −1 (Fig. (5).
The upper interface of the DSC is assumed to be planar although numerical simulations (Gerya et al., 2006) suggest that it may be corrugated by plumes developing in the DSC.A markedly non-planar interface would certainly produce scattering depending on length scale and shape.Our model of the DSC should be considered as an end-member model that avoids such complexity by starting from the simplest assumptions.

Effects of a deep subduction channel on wave propagation
Can a subduction channel presumably less than 10 km thick be detected by seismic methods?Tomographic studies appear not to be suitable because their resolution is generally insufficient to image such a thin channel.Methods using scattered or converted teleseismic waves (Bostock, 2002;Suckale et al., 2009) have potential to reach a resolution on the order of 10 km, but it is questionable whether they could distinguish the subduction channel from the oceanic crust.An alternative approach could be to search for specific phases in the seismic records which may serve as indicators for the existence of a channel.Since the channel is a small scale feature, visible effects on the seismic wave field can only be expected if the waves spend a considerable fraction of their propagation time inside the channel.For this reason, we consider here intermediate-depth earthquakes located in the subducted oceanic crust which directly radiate seismic energy into the subduction channel.Since both oceanic crust and channel are low-velocity regions for the major part of the propagation distance, we expect the formation of delayed high-amplitude guided waves.In the following, we will show by numerical simulations that a thin low-velocity channel on top of the subducted oceanic crust can significantly alter the seismic wave field observed at the surface.

Model setup
We use the 2-D Spectral Element Code by Komatitsch and Vilotte (1998) to simulate wave propagation in a 650 km wide and 300 km deep rectangular domain (Fig. 6).The domain is discretized into 835 × 96 quadrilateral elements with four internal Gauss-Legendre-Lobatto points per element.Source time function is a Ricker wavelet with a dominant source frequency of 2 Hz implying a significant spectral bandwidth from 1 to 5 Hz.We computed 120 000 time steps with a sampling interval of 0.0011 s.Synthetic seismograms are calculated for a line of 130 receivers at the surface spanning the entire width of the domain with a spacing of 5 km.
Receivers are numbered consecutively beginning with 1 from left (south) to right (north).Paraxial absorbing boundary conditions are applied at all boundaries except for the surface which is free.The computational domain has been chosen much deeper and wider than required to avoid disturbing reflections from the boundaries.The model is only 2-D implying excitation by line sources instead of point sources leading to reduced geometrical spreading and phase shifts.For comparison with observations, 3-D computations would be an option for the future.

Guided P waves
At first, we concentrate on the essential features of the compressional (P ) wave field.An explosive point source is put inside the oceanic crust at a depth of 108 km.Although an explosive source can not represent a real earthquake, it proves useful here because it does not radiate shear waves which, owing to their large amplitudes, hamper the proper visualization of details in the P wave field by means of snapshots.Moreover, an explosive source well approximates P waves from a double couple source with maximum radiation into the low-velocity waveguide formed by oceanic crust and DSC.This is confirmed by computations with a doublecouple source shown in context with S wave propagation below (Fig. 14).
The general propagation situation is displayed in Fig. 6.A P wave with nearly spherical wavefront has expanded from the source.In the oceanic crust, a guided wave has formed which is already delayed compared to the wavefronts in the surrounding mantle.In addition, some P S conversions have developed at the upper and lower boundary of the oceanic crust.
The wave field develops more structure while it passes through regions with significantly lower velocities than in the surrounding mantle.When the guided waves reach the bend of the slab and propagate into the serpentinized mantle wedge (Fig. 7) two distinct wave packages can be identified: a delayed high-amplitude guided wave train made up of the delayed direct P wave followed by PSP-conversions created at the crust-mantle interface, and an advancing re-Fig.6: Snapshot of vertical displacement illustrating general setup of the numerical model.Propagation time is 13.2 s.An explosive source radiates P waves travelling up the slab.Inside the oceanic crust, the P wave is already slightly delayed but the wavefront is still attached to the waves propagating in the surrounding mantle.At the upper and lower interface, PS conversions have developed.They can be identified by their smaller wavelength.Vertical axis annotation gives depth in km.Horizontal axis specifies index of receivers.Color shading for displacement displays negative values as blue and positive values as red shades.Color scale is nonlinear.RGB values are proportional to normalized absolute value of displacement to the power of 0.3 enhancing low-amplitude displacements.Displacements less than 1 percent of the maximum value are suppressed.
fracted P wave generated by the fast wave travelling in the high-velocity mantle below the interface.
Later, these wave trains appear as prominent arrivals at the surface.The models with and without DSC mostly differ with respect to the guided waves while the direct P wave (2) and the slab-refracted P wave (1) appear in both models.In the model without DSC, there is only one guided wave radiating into the wedge (3) and one attached to the oceanic crust (4), whereas in the model with DSC we find two guided waves emanating from the upper boundary of the DSC (3 and 3 ) and two guided waves attached to the lower boundary of the oceanic crust (4 and 4 ).
The snapshots shown in Fig. 8 allow a better understanding of the resulting signals at the surface receivers.In the first one, at 24.2 s after origin time (Fig. 8a), the guided waves reach the surface for the first time.The direct P wave is still the first arrival but about to be overtaken by the slab refracted P wave.In the model with DSC (left panel), the two guided wave trains 3 and 3 follow the direct wave with high amplitudes.Phases 4 and 4 are visible as weak signals in the oceanic crust.In the model without DSC (right panel), phase 3 generates a surface arrival while guided phase 4, attached to the oceanic crust, has developed a strong straight wavefront inside the wedge.Moreover, a secondary phase (5) has formed during travel through the serpentinized mantle wedge.
At 27.5 s after origin time, (Fig. 8b), the direct P wave (2) has just been overtaken by the slab refracted wave (1).In the model with DSC (left panel), the guided phases 3 and 3 still www.solid-earth.net/5/141/2014/Solid Earth, 5, 141-159, 2014 produce the largest signals while in the model without DSC (right panel) phase 4 becomes dominant and the secondary phase 5 starts to split off the main guided phase 4. At 36.3 s after origin time (Fig. 8c), the slab-refracted wave (1) has left behind the slower waves and forms the first arrival.The direct P wave (2) has been overtaken by the guided waves and is submerged by the strong phase 4. In the model without DSC (right panel), phase 4 now forms the major onset closely followed by the slower phase 5.A PS-conversion of phase 5 at the boundary between crust and serpentinized mantle wedge has also reached the surface as phase 7.In the model with DSC (left panel), the guided phases 3 and 3 are still dominant but partly overtaken by the guided phases 4 and 4 which were previously attached to the oceanic crust.The latter only form weak arrivals in the seismograms compared to the strong phase 4 in the model without DSC.The PS-conversion of phase 3 at the crustwedge boundary (phase 7) is still very weak.In addition, a converted PS-refraction from the slab starts to move towards the surface (phase 6) and a reflection at the slab interface of the surface-reflected direct P wave develops (phase 8).
With help of the snapshots, synthetic seismogram sections computed for the receivers at the surface (Fig. 9 and 10) can be readily interpreted.Remarkably, though the velocity models just differ by the presence or absence of the thin DSC, these record sections differ in the major arrivals.The pattern of onsets in the sections indeed define tell-tale signatures for the presence or absence of a DSC.The "no-DSC" signature (Fig. 9) is characterized by the dominating phase 4 originating as guided wave from the oceanic crust and its follow-up phase 5 that emerged during propagation through the slow mantle wedge.An arrival between phase 2 and 3 is missing.On the other hand, the "DSC"-signature is characterized by the strong phase pair 3 and 3 originating as guided waves from the DSC.Phase 4 is essentially missing and replaced by the weak phase pair 4 and 4 .It should also be noted that in the model with DSC the combination of direct wave (1), slab refracted (2) and guided waves (3, 3 ) forms uninterrupted wave trains of about 4 s length between receivers 60 and 70, a location corresponding to Crete in our subduction zone model.
To investigate whether the tell-tale signatures with respect to the presence or absence of a DSC are preserved when the source is moved up-or downwards inside the oceanic crust, we performed numerical simulations for varying source positions inside the oceanic crust (Fig. 11) for the model with DSC.If the source is too shallow, the phases 2, 3 and 3 merge and are no longer distinguishable.In our model, this occurs at source depths less than 70 km.For greater source depths, the "DSC"-signature can be clearly recognized.The additional phase 3 is present in all cases.The deeper the position of the source, the longer are the guided wave trains.This behavior is a consequence of the longer travel distance inside the waveguide formed by oceanic crust and DSC.For a source depth of 155 km the arrivals created by the guided waves become extremely long.Wave trains of up to 10 s length can be generated without any need for small-scale scatterers inside the medium.Due to the low velocity in the DSC, multiple P reflections and PS-conversions can build up in the waveguide and are able to conserve their amplitudes because of the velocity contrasts on both sides.Without DSC, there is no real difference between a deep or a shallow source with respect to the length of the guided wave.This behavior is caused by the steady velocity increase inside the oceanic crust with depth due to progressive transformation of basalts into less water-bearing, denser and faster mineral phases.The oceanic crust alone starts to become an efficient waveguide only at rather shallow depths.
What happens if the source is not inside the oceanic crust but either in the mantle below or in the DSC above?We performed two additional simulations where we moved the source from its original location within the oceanic crust at 108 km depth to a depth of 118 km in the mantle part of the slab and to a depth of 98 km in the DSC, respectively (Fig. 12).If the source is in the slab mantle, we observe a clear "no-DSC"-signature despite the presence of the DSC in the model.If the source is moved into the DSC, the record section is still of "DSC"-type.This finding may prove unfavourable for a detection of a DSC from seismic records.Earthquakes may well move into the mantle part of the slab because the temperature minimum also moves there during slab descent.On the other hand, if intermediate depth earthquakes are caused by dehydration reactions, it is highly probable that they occur in the oceanic crust which contains most of the subducted water and from where fluids are discharged into the overlying mantle wedge (Hacker et al., 2003a).
One could argue that errors in earthquake location may inhibit the identification of a DSC, for example, when a "no-DSC"-signature is observed for an earthquake erroneously located in the oceanic crust.Hence, seismic evidence for a DSC should always relate to observation of both signatures.If both are present, in the ideal case this could even be used to discriminate between earthquakes located within oceanic crust and DSC and those located deeper in the mantle part of the slab.
For efficiency reasons, the numerical simulations were carried out neglecting attenuation of seismic waves.To ensure the validity of the conclusions derived from the purely elastic simulations, we performed one simulation for a model with DSC assuming a spatial distribution of the quality factor as depicted in Fig. 4. The resulting record section for a source at 108 km depth does not significantly differ from the one shown in Fig. 10.In particular, the "DSC"-signature is fully preserved with attenuation.

Guided shear waves
Can shear waves help distinguish a subduction zone with DSC from one without?We performed further simulations for a double-couple source with a rupture plane striking     1) P wave refracted at the interface between subducted oceanic crust and lithospheric mantle, (2) direct P phase, (3) first guided P phase emerging from the deep oceanic crust, (4) second guided P phase emerging from the deep oceanic crust, (5) third guided P phase developing in the serpentinized mantle wedge, (6) refracted S phase converted from direct P phase at interface between subducted oceanic crust and lithospheric mantle, (7) S conversion of upgoing guided phase 5 at serpentinized wedge-crust interface, (8) surface reflection of direct P phase reflected back at slab interface.Solid lines crossing seismograms indicate lines of constant time for comparison with snapshots a-c of Fig. 8 on the right.Associated times from top to bottom are t = 24.2s, 27.5 s and 36.3 s.Source is explosive and at 108 km depth.perpendicular to the subduction direction (profile direction) and dipping 15 degrees downwards.Since the deeper part of the slab dips at 30 degrees, the angle of 15 degrees was chosen to achieve both sufficient P and S wave radiation into the oceanic crust.The source depth is again 108 km as in our reference case for P waves.The synthetic record sections are slightly more complicated because SP conversions arrive now before the guided shear waves.As for the P waves, again tell-tale signatures can be observed in the record sections which allow a discrimination between models with and without DSC by one quick glance (Fig. 13).The guided shear wave signals appear within the rectangular boxes inserted into the record sections.Beyond receiver number 65, three distinct phases (3a, 3b, 3c) can be recognized for the model with DSC whereas there is only one for the model without.At receivers 40 to about 65, the situation reverses: the model without DSC now exhibits a package of three phases while there is only one in the model with DSC.The phases preceding those in the rectangular box appear for both models Numbers from 1 to 7 mark major seismic phases.(1) P-wave refracted at the interface between subducted oceanic crust and lithospheric mantle, (2) direct P-phase, (3) and (3') two guided P-phases emerging from the deep oceanic crust and radiated into the serpentinized mantle wedge, ( 4) and (4') guided P-phases emerging from the deep oceanic crust staying attached to the oceanic crust, (6) refracted Sphase converted from direct P-phase at interface between subducted oceanic crust and lithospheric mantle, ( 7  and are not suitable for discrimination purposes.They are the direct S wave (1), a slab-refracted S wave (2), a slabrefracted SP surface conversion (4) and a slab-reflected SP conversion (5).The snapshots in Fig. 14 provide some additional understanding of guided shear wave propagation: in the model with DSC, an extended M-shaped guided wave package develops in the steeper part of the slab from which phases 3a-c develop later.As for P waves, the guided wave train in the model without DSC is rather compact (Fig. 14a).Once these waves move into the serpentinized mantle wedge, they develop strong wavefronts (Fig. 14b) which keep their identity while passing through the wedge and arrive as strong individual phases at the surface (Fig. 14c).In the former case, the record section exhibits a "no DSC"-signature (compare to Fig. 9), in the latter case, it exhibits a clear "DSC" signature (compare to Fig. 10).Phases are numbered consistent with Figs. 9 and 10.

Observations of guided
EGELADOS experiment (Friederich and Meier, 2008) in the southern Aegean.The seismograms are band-passed filtered in the frequency range from 2 to 5 Hz in order to mimic the effect of the Ricker wavelet used in the numerical modeling.We selected seismograms from stations located on Crete for the most part, complemented by some low-noise oceanbottom hydrophone records.The data records are sorted ac-cording to epicentral distance.The P signals are characterized by multiple high-amplitude wave trains of several seconds length that are also found in the synthetic seismograms where they are produced by guided waves in oceanic crust and DSC.
As the installed stations are not aligned along a section plane containing the hypocenter, and their spacing is too large and irregular, the available data are not suitable for recognizing the DSC or no-DSC signature identified within our simulations.New dedicated seismic field experiments are required to clarify this issue.The MEDUSA experiment (Suckale et al., 2009) carried out in the Ionian region may provide the required station density and alignment to reveal DSC signatures.
Long wave trains could alternatively be generated by scatterers of a size on the order of a wavelength.However, numerical simulations done by the authors (unpublished) indicate that the length and amplitude of the observed P wave trains can only be reproduced by placing strong, kilometersized velocity anomalies of 20 to 40 percent contrast into the oceanic crust.The beauty of our model is that the long, high amplitude wave trains can be effectively generated by a simple, thin, homogeneous low-velocity channel on top of the oceanic crust.

Discussion
The fundamental question is to which extent the results of our numerical simulations can be applied to natural situations.More specifically, can seismic records that display the "DSC signature" be taken to indicate the presence of a subduction channel, or, vice versa, can the absence of this signature in seismic records be taken to discard the hypothesis?At this point, it must be emphasized that the structural model of a DSC underlying our simulations is extremely simplistic, although it incorporates presumably realistic structural elements, phase assemblages and elastic properties based on mineral physics.These properties are poorly constrained, mainly due to uncertainties with respect to bulk composition and mineral assemblage, widespread disequilibrium and metastability, heterogeneity on all length scales, probably steep temperature gradients and possibly complex internal structures which can only tentatively be inferred from geological observations in exhumed high P /T metamorphic terranes, leaving ample space for interpretation.
Moreover, at present it seems impossible to identify classes or styles of DSCs, as a "quenched" DSC structure is not available for analysis at the surface.Numerical simulations are probably the only way to depict possible histories of DSCs, choosing starting and boundary conditions which at least appear realistic, corresponding to common sense.This situation could improve when presently active subduction channels at depth can be identified by seismological means.Only in this case, true snapshots may be correlated with current state and history of the respective subduction system, and eventually open the chance to identify different types of DSC.Our present contribution is thought to provide a small step into this direction.
Limited by the wavelength of recordable seismic signals and the resolution of seismic experiments, our model of a subduction channel accordingly integrates over a large volume.Most importantly, it substitutes highly heterogeneous material in a real subduction channel by a homogenized, effective material which neglects heterogeneities below the kilometer-scale and temperature gradients within the subduction channel.The temperature field is expected to be con-trolled by the overall flow pattern which may change in space and time with a potentially strong feedback.In addition, our model adopts a simple crustal structure and is only twodimensional.In contrast, the Hellenic subduction zone (HSZ) to which it is applied exhibits considerable curvature and a possibly corrugated slab shape (Knapmeyer, 1999;Papazachos et al., 2000).Finally, anisotropy inside the subduction channel is neglected, as many exposed high P /T metamorphic rock assemblages show structural patterns which cannot necessarily be expected to produce anisotropy on the length scale resolved by seismic experiments.The same holds for BIM structures taken to represent a subduction channel here.
Available data from the HSZ are not suitable to search for "DSC signatures" by tracking seismic phases over epicentral distance, because the seismic stations in the forearc of the HSZ are not located on a straight profile but exhibit varying azimuths with respect to the intermediate-depth hypocenters (map in Fig. 15d slab, reaching from the outer rise to the backarc, could remedy the situation.In such an experiment, stations should be linearly aligned with a spacing of less than 5 km.Even though intermediate-depth earthquake data generally exhibit more complicated waveforms than those obtained from corresponding simulations, the characteristic guided P wave trains of several seconds length are well discernible.These strong signals, tracked and correlated over epicentral distance, would allow identification of tell-tale signatures for the presence or absence of a subduction channel.In the Aegean, as probably in most subduction zones worldwide, ocean bottom seismometers would be required besides land stations to provide the required data. The HSZ is considered an attractive target region for this type of experiment, as an active DSC has been proposed based on independent evidence (Meier et al., 2007).These authors found that the outline of a region lacking microseismicity at a depth of about 20 to 40 km coincides with the coast line of Crete.Significant uplift of Crete since 4 Ma in an overall extensional environment is therefore suspected to be correlated with anomalous material at depth.Based on position within the HSZ, the properties are tentatively proposed to reflect extrusion of material from a DSC, which became fully established at about 4 Ma after Miocene microcontinent collision (Meier et al., 2007).The fact that only discrete parts of the forearc, forming the islands, undergo uplift may be explained by disintegration of the DSC into distinct fingers, as the radius of the curved plate interface increases towards the surface.If signatures of a DSC could be observed in a seismic experiment at the HSZ, they could be taken to add another piece of evidence in favour of the above hypothesis..A limiting factor for any such passive seismic experiment is the occurence of suitable intermediate depth earthquakes.For example, observations of seismicity in the Aegean (Hatzfeld and Martin, 1992) indicate that intermediate depth earthquakes cluster in a small region close to the islands Nisyros and Astypalea in the southeastern corner of the Aegean.Using these earthquakes would imply a profile across Karpathos, for which the layout requires a large share of ocean bottom seismometers.For a more convenient seismic profile across Crete, intermediate depth earthquakes further to the west along the strike of the subduction zone would be prerequisite, which are rarely observed.As a compromise, a profile oriented towards the Nisyros earthquake cluster could be feasible (green line in Fig. 15d), albeit being oblique to the dip of the slab.In that case, one would have to cope with three-dimensional propagation effects owing to slab curvature.

Conclusions
We performed numerical simulations of seismic wave propagation in a generic subduction zone model to explore potential seismological evidence for an active deep subduction channel.The distribution of seismic velocities in our model is based on geometric and structural information, predicted temperature field, phase diagrams for MORB and harzburgite, and properties of minerals at high pressures and temperatures predicted using equations of state.
Although suffering from various uncertainties, the model produces self-consistent values for P , S velocity and density and their variation with pressure and temperature.It produces velocity gradients in the slab and the mantle wedge and a steady velocity increase in the oceanic crust, an essential feature for the development of guided waves.It thus represents a more realistic model of a subduction zone than models assuming subregions with constant velocities gleaned from tomographic images.
The deep subduction channel itself is modeled as a thin layer located on top of the subducted oceanic crust.This layer is composed of eclogitic blocks floating in a serpentinite matrix.Homogenization calculations show that such a BIM structure exhibits frequency-dependent seismic velocities with values intermediate between those of block material and pure matrix.Since seismic velocities of serpentinite are low compared to those of typical mantle rock, and the volume fraction of eclogitic blocks is small (between 5 and 50 %), the subduction channel is predicted to represent a low-velocity zone acting as wave guide.Petrophysical calculations suggest that seismic velocities of serpentinite matrix are nearly constant over the entire length of the subduction channel, the effects of temperature and pressure approximately cancelling each other.
Our numerical simulations focus on high-frequency seismic waves emitted from intermediate depth earthquakes.Analysis of synthetic seismograms for receivers located in the forearc indicate significant differences of P and S waveforms for models with and without deep subduction channel.Although only a very small volume of the model is covered by the subduction channel, the presence or absence of a channel is found to affect major seismic arrivals.Record sections along profiles crossing the forearc parallel to slab dip exhibit tell-tale patterns which allow to discriminate between models with and without deep subduction channel.The "DSC" and "no-DSC" signatures are found to represent stable features which remain preserved when the earthquake is shifted up or down within the oceanic crust.The signature also remains visible when the source mechanism is changed.Only in case the earthquake occurs below the oceanic crust, the "DSC" signature is lost and a "no-DSC" signature appears.Our results suggest that the hypothesis of deep subduction channel, explaining rapid exhumation of high P /T metamorphic rocks, might be tested by an appropriately designed seismic experiment.
Ideally, such an experiment would provide observations of both DSC and no-DSC signatures from different earthquakes.This way, it would become unnecessary to discriminate between earthquakes located within DSC or oceanic crust and those located in the slab mantle, but the position could be inferred from the signature.

Fig. 2 :
Fig. 2: (Top) Regions within which different laws for calculating temperature were applied, overlain by black lines indicating structural regions of subduction zone model derived from seismological studies, (bottom) temperature distribution on a vertical section through subduction zone model.Structural regions now indicated by white lines.Temperature at the base of the lithosphere has been assumed as 1250 • C. Units: degrees Celsius

Fig. 2 :
Fig. 2: (Top) Regions within which different laws for calculating temperature were applied, overlain by black lines indicating structural regions of subduction zone model derived from seismological studies, (bottom) temperature distribution on a vertical section through subduction zone model.Structural regions now indicated by white lines.Temperature at the base of the lithosphere has been assumed as 1250 • C. Units: degrees Celsius.

Fig. 2 .
Fig.2.At the grid nodes, lithostatic pressure is calculated using the vertical density distribution of a standard earth model.Temperature at the grid nodes is calculated as described above (see Fig.2).Depending on the bulk compositions at a given grid point (continental crust, MORB, harzburgite, dry harzburgite), we look up the stable phase and its modal composition in the corresponding phase diagram (Fig.3) and calculate elastic moduli and density for each constituting mineral using the theory of SLB2005 and the corresponding parameters in the mineral database.Effective elastic moduli of the phase are obtained by a volume fraction weighted average of the individual reciprocal elastic moduli (Reuss average).Density is calculated by averaging the individual density values according to the volume fractions.This procedure is certainly not satisfactory as it is known that the effective elastic moduli may have varying values depending on the spatial distribution of the modal constituents inside the phase.We have chosen to display here the Reuss averaged values to give an impression of the variations of seismic wave speeds inside a subduction zone.Hashin-Shtrikman bounds(Hashin and Shtrikman, 1963) can be readily calculated and provide an estimate of the variability of these values.Since, however, seismic wave speeds derived from tomographic imaging represent averages over rock volumes of kilometer size, a homogenization scheme is needed to obtain accurate estimates for the effective elastic moduli of such a rock volume.This is particularly true in subduction zones where large changes of temperature and also stress may occur within small distances.Clues of the spatial distribution of rocks inside a seismically

Fig. 4 :
Fig.4: Distribution of P-wave velocities (top left) and S-velocities (top right) in km/s, density in g/cm 3 (bottom left) and quality factor Q (bottom right), calculated from the temperature model using SLB2005 theory, the phase diagrams ofHacker (2003)  and our mineral data base which was derived from those ofHolland andPowell (1998) and Hacker (2003).

Fig. 5 :
Fig. 5: Distribution of P-wave velocity and S-wave velocity for the subduction zone model including a DSC

Fig. 4 :
Fig.4: Distribution of P-wave velocities (top left) and S-velocities (top right) in km/s, density in g/cm 3 (bottom left) and quality factor Q (bottom right), calculated from the temperature model using SLB2005 theory, the phase diagrams ofHacker (2003)  and our mineral data base which was derived from those ofHolland andPowell (1998) and Hacker (2003).

Fig. 5 :
Fig. 5: Distribution of P-wave velocity and S-wave velocity for the subduction zone model including a DSC Fig. 5: Distribution of P wave velocity and S wave velocity for the subduction zone model including a DSC.

Fig. 7 :
Fig. 7: Snapshot of P wave propagation for a model with DSC (left) and without DSC (right).Propagation time is 18.15 s.Phase 1 is the slab refracted wave, phase 2 the direct P phase and 3 and 4 are guided waves.Note that in the model with DSC there are two guided waves (3, 3 and 4, 4').Phases 4 and 4 in the model with DSC are much weaker and phases 3 and 3 are much stronger than their counterparts in the model without DSC.Numbering of phases is the same as in Fig. 9. See Fig. 6 for explanation of axis labels and color shading.

Fig. 8 :
Fig. 8: Sequence of snapshots illustrating the evolution of the wave field.Left column: model with DSC; Right column: model without DSC.Vertical axis annotation gives depth in km.Horizontal axis specifies index of receivers.(a) t = 24.2s, guided phases 3 have reached the surface.The direct phase still forms the first onset.(b) t = 27.5 s, slab-refracted wave 1 overtakes direct wave 2; phase 4 forms dominant signal in model without DSC and a secondary guided phase 5 develops in the slow serpentinized wedge; on the contrary, in the model with DSC, phases 3 form distributed signals at the surface and phases 4 stay weak.(c) t = 36.3s.In the model with DSC, phases 4 overtake phases 3 and form weak precursory signals; an S conversion of phase 5 reaches the surface (phase 7).The slab refracted S conversion (phase 6) becomes visible and a reflection from the slab interface of the P surface reflection emerges.Strongest signal in model with DSC is phase 3. Phase 4 is strongest signal in model without DSC.Numbering of phases is the same as in Fig. 9. See Fig. 6 for explanation of axis labels and color shading.

Fig. 9 :
Fig. 9: Model without DSC: time-reduced record section of synthetic horizontal component seismograms computed at receivers at the surface.Numbers from 1 to 8 mark major seismic phases.(1)P-wave refracted at the interface between subducted oceanic crust and lithospheric mantle, (2) direct P-phase, (3) first guided P-phase emerging from the deep oceanic crust, (4) second guided P-phase emerging from the deep oceanic crust, (5) third guided P-phase developing in the serpentinized mantle wedge, (6) refracted S-phase converted from direct P-phase at interface between subducted oceanic crust and lithospheric mantle, (7) S-conversion of upgoing guided phase 5 at serpentinized wedge-crust interface, (8) surface reflection of direct P phase reflected back at slab interface.Solid lines crossing seismograms indicate lines of constant time for comparison with snapshots a-c of Fig.8on the right.Associated times from top to bottom are: t = 24.2s, 27.5 s and 36.3 s.Source is explosive and at 108 km depth.

Fig. 9 :
Fig. 9: Model without DSC: time-reduced record section of synthetic horizontal component seismograms computed at at the surface.Numbers from 1 to 8 mark major seismic phases.(1) P wave refracted at the interface between subducted oceanic crust and lithospheric mantle, (2) direct P phase, (3) first guided P phase emerging from the deep oceanic crust, (4) second guided P phase emerging from the deep oceanic crust, (5) third guided P phase developing in the serpentinized mantle wedge, (6) refracted S phase converted from direct P phase at interface between subducted oceanic crust and lithospheric mantle, (7) S conversion of upgoing guided phase 5 at serpentinized wedge-crust interface, (8) surface reflection of direct P phase reflected back at slab interface.Solid lines crossing seismograms indicate lines of constant time for comparison with snapshots a-c of Fig. 8 on the right.Associated times from top to bottom are t = 24.2s, 27.5 s and 36.3 s.Source is explosive and at 108 km depth.

Fig. 10 :
Fig. 10: Model with DSC: time-reduced record section of synthetic horizontal component seismograms computed at receivers at the surface.
) S-conversion of upgoing guided phases 3 at serpentinized wedge-crust interface.Solid lines crossing seismograms indicate lines of constant time for comparison with snapshots a-c of Fig. 8 on the left.Associated times from top to bottom are: t = 24.2s, 27.5 s, 36.3 s.Source is explosive at 108 km depth.

Fig. 10 :
Fig. 10: Model with DSC: time-reduced record section of synthetic horizontal component seismograms computed at receivers at the surface.Numbers from 1 to 7 mark major seismic phases.(1) P wave refracted at the interface between subducted oceanic crust and lithospheric mantle, (2) direct P phase, (3) and (3') two guided P phases emerging from the deep oceanic crust and radiated into the serpentinized mantle wedge, (4) and (4') guided P phases emerging from the deep oceanic crust staying attached to the oceanic crust, (6) refracted S phase converted from direct P phase at interface between subducted oceanic crust and lithospheric mantle, (7) S conversion of upgoing guided phases 3 at serpentinized wedge-crust interface.Solid lines crossing seismograms indicate lines of constant time for comparison with snapshots a-c of Fig. 8 on the left.Associated times from top to bottom are t = 24.2s, 27.5 s, 36.3 s.Source is explosive at 108 km depth.
P waves in the Aegean P phases of typical seismograms from intermediate-depth earthquakes recorded in the forearc of the Hellenic subduction zone are shown in Fig. 15.The data were acquired between October 2005 and March 2007 during the

Fig. 11 :
Fig. 11: Time-reduced synthetic record section for the model with DSC for different source depths inside the oceanic crust.Depths: (a) 88 km, (b) 108 km, (c) 128 km, (d) 155 km.The additional phase 3 is present for all depths.Note the increasing length of the guided wave train with increasing source depth.

Fig. 12 :
Fig.12: Time-reduced synthetic horizontal component record sections for the model with DSC and explosive sources located below the oceanic crust (left panel) and inside the DSC (right panel).In the former case, the record section exhibits a "no DSC"-signature (compare to Fig.9), in the latter case, it exhibits a clear "DSC" signature (compare to Fig.10).Phases are numbered consistent with Figs. 9 and 10.

Fig. 13 :
Fig. 13: Time-reduced synthetic record sections for a double-couple source at 108 km depth showing S waves.Left panel: model with DSC.Right panel: model without DSC.Marked phases: (1) direct S wave, (2) S wave refracted at slab interface, (3) guided wave, (3a-c) multiple guided waves in model with DSC, (4) SP surface conversion refracted at slab interface, (5) SP surface conversion reflected at slab interface.Phases 1, 2, 4 and 5 appear in models with and without DSC.Note the three distinct guided waves in the model with DSC.Boxes enclose signals related to guided waves.Lower panels: zoom into the receiver range 62 to 75 where the three guided arrivals are best visible.

Fig. 14 :
Fig. 14: Sequence of snapshots of vertical displacement for double couple source.Left panels: model with DSC.Right panels: model without DSC.Vertical axis annotation gives depth in km.Horizontal axis specifies index of receivers.Phase identification is given in Fig. 13.(a) t = 24.75s.Guided S waves about to enter the serpentinized mantle wedge; several wave packets have developed in model with DSC (phases 3a, 3b, 3c).(b) t = 34.1 s.Guided wave already in serpentinized area, direct S wave already at surface, note strong SP surface conversion.(c) t = 44 s.Guided waves about to reach the surface.See Fig. 6 for explanation of axis labels and color shading.

Fig. 15 :
Fig. 15: Examples of time-reduced seismograms for intermediate depth events recorded at stations located in the forearc of the Hellenic subduction zone as shown in the event and station map in (d).(a) Event e0565.217.06,depth 87 km; (b) event e1049.334.05,depth 99 km; (c) event e1866.350-06,depth 110 km.Note the multi-onset character of the P phases and the generally long complex wave trains.Red diamonds in (d) show event epicenters and green line indicates possible profile for a seismic experiment.

Table 1 :
List of constants used to compute temperature and pressure dependent values for the inverse quality factor according to Eq. 1 Bezacier et al. (2010)s 6.0 km s −1 , 2.73 km s −1 and 2680 kg/m 3 , respectively.Recently,Bezacier et al. (2010)obtained isotropic Reuss average aggregate P and S velocities of antigorite using Brillouin-scattering of 6.17 km s −1 and 3.32 km s −1 .While the value for P velocity is consistent with Hacker's value, S velocity ofBezacier et al. (