Chapter 5: Geodynamics

Chapter 5 Geodynamics 5.1 Heat flow Thermally controlled processes within Earth include volcanism, intrusion of igneous rocks, metamorphism, convection within the mantle and outer core, and plate tectonics. The global heat flow can be measured by measuring the temperature gradient everywhere at the surface of the Earth. This gives us an estimate of the mean rate of heat loss of the Earth, which can be broken up into various components (Table 7.3 in Fowler): Continents Oceans Area (km2 ) 201 309 Conductive cooling Hydrothermal circulation Total Earth 510 Heat Flow (mWm−2 ) 58 100 [66] [34] 83 Heat Loss (1012 W) 11.5 30.4 [20.3] [10.1] 41.9 The amount of heat lost through the ocean basins in enormous! — up to 73%! (The oceans cover about 60% of the Earth’s surface). This was a famous paradox before the discovery of plate tectonics. It was well known that the abundance of radioactive elements (which are a source of heat through radioactive decay) in the ocean basins was much lower than that in the continents. So what causes the significantly higher heat flow in the oceans? With the discovery of plate tectonics it was realized that most of the heat loss occurs through the cooling and creation of oceanic lithosphere. The mean rate of plate generation therefore depends on the balance between the rate of heat production within the Earth and the rate of heat loss at the surface. In this course we will address some of the basic concepts of heat flow and Earth’s thermal structure, and we will discuss in some detail the cooling of oceanic lithosphere and the implications of Earth thermal structure for mantle convection. Heat sources There are several possibilities for the source of heat within the earth: 1. ”Original” or ”primordial” heat; this is the release of heat due to the cooling of the Earth. The amount of heat released by this process can be estimated by calculating the heat released by a change in temperature of 1◦ at constant pressure. This depends on the specific heat, CP which is the energy that is needed to heat up 1 kg of material by 1◦ (i.e., it’s a material property). We can do a quick calculation to find out how much heat would be released by dropping the temperature of the mantle by 1◦ C (Let’s for now ignore latent heat due to phase changes): • Mantle; for silicates:CP = 7.1 × 102 Jkg−1◦ C−1 ; the mass of the mantle is about 4.1×1024 kg • Core; for iron: CP = 4.6 × 102 Jkg−1◦ C−1 ; the mass of the core is about 1.9×1024 kg For ΔT = 1◦ C this gives ΔE = 3.7×1027J. In absence of any other sources for heat production, the observed global heat flux of 4.2×1013W can thus be maintained by a cooling rate of 4.2×1013 [W] divided by 3.7×1027 [J]= 1.1×10−14◦ Cs−1 . In other words, since the formation of Earth, 4.5 Ga ago, the average temperature would have dropped by ΔT ≈ 1, 500◦C. Note that the actual cooling rate is much lower because there are sources of heat production. 2. Gravitational potential energy released by the transfer of material from the surface to depths. Imagine dropping a small volume of rock from the crust to the core. The gravitational potential energy released would be: ΔE = Δρgh, with g ≈ 10ms−2 and h = 3 × 106 m ρsilicates ≈ 3 × 103 kgm−3 and ρiron ≈ 7 × 103 kgm−3 , so that Δρ = 4 × 103 kgm−3 . ΔE ≈ 1.2 × 1011 Jm−3 . The present-day heat flux would thus be equivalent to dropping a volume of about 350 m3 every second. This is equal to dropping a 22 m thick surface layer every million years. So even if a small amount of net differentiation were taking place within the earth, this would be a significant source of heat! 3. Radioactive decay: for an order of magnitude calculation, see Stacey 6.3.1. The bottom line is that for the Earth a very significant fraction of heat loss can be attributed to radioactive decay (primarily of Uranium (U), Thorium (Th) and Potassium (K). More, in fact, than can be accounted for by heat production of the MORB source. Heat transfer The actual cooling rate of the Earth depends not only on these sources of heat, but also on the efficiency at which heat is transferred to and lost at the Earth’s surface. How does heat get out of the system? Conduction — this will be discussed below in the context of the cooling of oceanic lithosphere. Convection — For example, in the mantle and core. Radiation — most of the heat that the Earth receives from external sources (i.e. the Sun) is radiated out. Radiation The net effect is that the Earth is cooling at a small rate (of the order of 50-100◦C per Ga!) (See Stacey (1993), p. 286.) Figure 5.1: 5.2 Heat flow, geothermal gradient, diffusion The rate of heat flow by conduction across a thin layer depends on 1. the temperature contrast across the layer (ΔT ) 2. the thickness of the layer (Δz) 3. the ease with which heat transfer takes place (which is determined by the thermal conductivity k). The thinner the layer and the larger the temperature contrast (i.e., the larger the gradient in temperature), the larger the heat flow. In other words, the heat flow q at a point is proportional to the temperature gradient at that point. This is summarized in Fourier’s Law of conduction: q = −k∇T ≈ −k ΔT z ˆ Δz (5.1) where the minus sign indicates that the direction of heat flow is from high to low tempertaures (i.e., in the opposite dirtection of z if z is depth.). (For simplicity we talk here about a 1D flow of heat, but Fourier’s Law is also true for a general 3D medium). We can use this definition to formulate the conduction (or diffusion) equation, which basically describes how the temperature per unit volume of material changes with time. This change depends on 1. the amount of heat that flows in or out of the system which is described by the divergence of heat flow 2. the amount of heat produced within the volume (denoted by the density of heat sources A) 3. the coupling between this change in heat and a change in temperature (which is controlled by the specific heat) The thermal diffusion equation is given by: ρCP ∂T = −∇ · q + A ∂t (5.2) Or: the change in heat content with time equals the divergence of the heat flow (into and out of the volume) and the generation of heat within the volume. Combined with Fourier’s Law the diffusion equation can be written as ρCP ∂T = −∇ · (−k∇T ) + A = k∇2 T + A ∂t (5.3) In a situation of steady-state the diffusion equation transforms to the expression of the geotherm, the variation of temperature with depth in the Earth: k∇2 T + A = ρCP ∂T A = 0 ⇒ ∇2 T = − ∂t k (5.4) If there is no heat production (by radioactive decay), i.e., A = 0, then the temperature increases linearly with increasing depth. If A �= 0 then the temperature/depth profile is given by a second-order polymomial in z. In other words, the curvature of the temperature-depth profile depends on the amount of heat production (and the conductivity). Figure 5.2: Heat production causes nonlinear geotherms. A typical value for the geotherm is of order 20 Kkm−1 , and with a value for the conductivity k = 3.0 Wm−1 K−1 this gives a heat flow per unit area of about 60 mWm−2 (which is close to the global average, see table above). If the temperature increases according to this gradient, at a depth of about 60 km a temperature of about 1500 K is reached, which is close to or higher than the melting temperature of most rocks. However, we know from the propagation of shear waves that the Earth’s mantle behaves as a solid on short time scales (µ > 0). So what is going on here? Actually, there are two things that are important: 1. At some depth the geothermal gradient is no longer controlled by conductive cooling and adiabatic compression takes over. The temperature gradient for adiabatic compression (i.e., the change of temperature due to a change of pressure alone, without exchange of heat with its invironment) is much smaller than the gradient in the conducting thermal boundary layer. 2. With increasing pressure the temperature required for melting also increases. In fact it can be shown that with increasing depth in Earth’s mantle, the actual temperature increases (from about 0◦ C at the surface to about 3,500 ± 1000◦ C at the core-mantle boundary CMB) but the melting temperature Tm increases even more as a result of the increasing pressure. Consequently, at increasing depth in the mantle the ratio of T over Tm (the homologous temperature) decreases. At even larger depth, in Earth’s core, the temperature continues to increase, but the melting temperature for pure iron drops (pure chemical compounds — such as pure iron — typically have a lower melting temperature then most mixtures — such as silicate rock) so that the actual temperaute ex ceeds the melting temperature and the material is in liquid state. Even tough the mantle is ’solid’ it behaves as highly viscous fluid so that flowis possible over very long periods of time. 5000 5000 4000 Solidus 3000 3000 Temperature 2000 2000 Temperature (oC) Temperature (K) 4000 1000 1000 0 0 0 A 2000 D'' 6000 4000 Depth (km) 6371 5150 2891 400 670 L Inner Core Outer Core Lower Mantle L = Lithosphere (0-80 km) A = Asthenosphere (80-220 km) D'' = Lower-Mantle D'' Layer 400, 670 = Phase Transitions Figure by MIT OCW. Figure 5.3: Geotherms in the Earth. If we ignore heat production by radioactive decay we can simplify the conduction equation to ρCP k ∂T ∂T ∇2 T = κ∇2 T = k∇2 T ⇒ = ρCP ∂t ∂t (5.5) with κ the thermal diffusivity κ= k ρCP (5.6) We will look at solutions of the diffusion equation when we discuss the cooling of oceanic lithosphere after its formation at the mid oceanic ridge. Before we do that let’s look at an important aspect of the diffusion equation. From a dimensional analysis of the diffusion equation ∂T = κ∇2 T ∂t (5.7) we see that the diffusivity κ has the dimension of length2 × time−1 . We can √ now define a diffusion length L as L = κt. If a temperature change occurs at some time t0 , then after√a characteristic time interval τ it will have ’propagated’ over a distance L = κτ through the medium with diffusivity κ. Similarly, it takes a time l2 /κ for a temperature change to propagate over a distance l. 5.3 Thermal structure of the oceanic lithosphere Introduction The thermal structure of the oceanic lithosphere can be constrained by the observations of: 1. Heat flow 2. Topography (depth of the ocean basins) 3. Gravity (density depends inversely on temperature) 4. Seismic velocities (µ = µ(T ), λ = λ(T )); in particular, surface waves are sensitive to radial variations in wave speed and surface wave dispersion is one of the classical methods to constrain the structure of oceanic (and continental) lithosphere. In the following we address how the heat flow and the depth of ocean basins is related to the cooling of oceanic lithosphere. The conductive cooling of oceanic lithosphere when it spreads away from the mid-oceanic ridge can be described by the diffusion equation ∂T = κ∇2 T + A ∂t (5.8) We will simplify this equation by (1) ignoring the heat production by radiocative decay, so that A = 0 (this is reasonable for the oceanic lithosphere since the basalts do not contain a significant fraction of major radio-isotopes Uranium, Potassium, and Thorium)1 , and (2) by assuming a 2D geometry so that we can ignore the variations in the y direction. The latter assumption is justified for regions away from fracture zones. With these simplifications the diffusion equation would reduce to a � � 2 ∂ T ∂2T ∂T (5.9) = κ∇2 T = κ + ∂t ∂x2 ∂z 2 with z the depth below the surface and x the distance from the ridge. The variation in temperature in a direction perpendicular to the ridge (i.e., in the spreading direction x) is usually much smaller than the vertical gradient. In that case, the heat conduction in the x direction can be ignored, and the cooling of a piece of lithosphere that moves along with the plate, away from the ridge, can be described by a 1D diffusion equation: ∂T =κ ∂t � ∂ 2T ∂z 2 � (5.10) (i.e., the ’observer’, or the frame of reference, moves with the plate velocity u = x/t). Note that, in this formulation, the time plays a dual role: it is used as the time at which we describe the temperature at some depth z, but this also relates to the age of the ocean floor, and thus to the distance x = ut from the ridge axis). u Ridge x=0 SURFACE u T = Ts x LITHOSPHERE ISOTHERM Tm ASTHENOSPHERE y Figure by MIT OCW. Figure 5.4: The cooling of oceanic lithosphere. Ridge T = Ts u T = Tm u q q t=0 x t= u u q t = t1 t = t2 Figure by MIT OCW. Figure 5.5: Bathymetry changes with depth. The assumption that the oceanic lithosphere cools by conduction alone is pretty good, except at small distances from the ridge where hydrothermal circulation (convection!) is significant. We will come back to this when we discuss heat flow. There is a still ongoing debate as to the success of the simple cooling model described below for large distances from the ridge (or, equavalently, for large times since spreading began). This is important since it relates to the scale of mantle convection; can the cooling oceanic lithosphere be considered as the Thermal Boundary Layer (across which heat transfer occurs primarily by conduction) of a large scale convection cell or is small scale convection required to explain some of the observations discussed below? See the recent Nature paper by Stein and Stein, Nature 359, 123–129, 1992. Cooling of oceanic lithosphere: the half-space model The variation of temperature with time and depth can be obtained from solving the instant cooling problem: material at a certain temperature Tm (or T0 in Turcotte and Schubert) is instantly brought to the surface temperature where it is exposed to surface temperature Ts (see cartoons below; for a full derivation, see Turcotte and Schubert). y=0 T0 T y=0 t = 0­ y Ts T0 T y=0 T0 t = 0+ Ts T t>0 y y Figure by MIT OCW. Figure 5.6: The heating of a halfspace Diffusion, or relaxation to some reference state, is described by error functions2 , and the solution to the 1D diffusion equation (that satisfies the appropriate boundary conditions) is given by � z T (z, t) = Tz (t) = Ts + (Tm − Ts ) erf √ 2 κt � (5.11) or � z T (z, t) − Ts = (Tm − Ts ) erf √ 2 κt � (5.12) with T (z, t) the temperature within the cooling boundary layer, Ts and Tm the temperature at the surface and in the mantle, respectively, κ the thermal 2 So called because they are integrations of the standard normal distribution. diffusivity3 , κ = k/ρCp (k is the thermal conductivity and Cp the specific heat), and the error function operating on some argument η defined as 2 erf(η) = √ π �η 2 e−u du (5.13) 0 The so called complementary error function, erfc, is defined simply as erfc(η) = 1 − erf(η). The values of the error function (or its complement) are often presented in table form4 . Figure 5.7 depicts the behavior of the error function: when the argument increases the function value ’creeps’ asymptotically to a value erf = 1. 1.0 0.8 erf � 0.6 0.4 erfc � 0.2 0 0 1 2 3 � Figure by MIT OCW. Figure 5.7: Error function and complimentary error function. Let’s look at the temperature according to (5.12) for different boundary conditions. For large values of z the solution of the diffusion equation becomes T (∞, t) = Tm ; at the surface, z = 0 so that T (0, t) = Ts , and= after a very long 3 The thermal diffusivity κ has the dimension of distance2 /time; a typical value for κ is 1 √ mm2 /s. The square root of the product κt is porportional to the diffusion length L ∼ κt. If the temperaure changes occur over a characteristic time interval t they will propagate a distance of the order of L. Similarly, a time l2 /κ is required for temperature changes to propagate distance l. 4 Type help erfin MatlabTM time, T (z, ∞) = Ts , i.e., the whole system has cooled so that the temperature is the same as the surface temperature everywhere. √ For the Earth we can set Ts = 0◦ C so that T ≈ Tm erf(η), η = z/(2 κt) for most practical purposes; but the above formulas are readily applicable to other boundary layer problems (for instance to the cooling of the lithosphere on Venus where the surface temperature is much than that at Earth). Examples of the geothermal gradient as a function of lithospheric age are given in the diagram below (from Davies & Richards, ”Mantle Convection” J. Geol., 100, 1992). Temperature (oC) 0 0 1000 5 Myr Depth (km) 25 Myr 100 Myr Dry solidus 100 Craton 200 Figure by MIT OCW. Figure 5.8: Cratonic and oceanic geotherms. Figure 5.9 (from Turcotte & Schubert, 1982) shows a series of isotherms (lines of constant temperature (i.e, T (z, t) − Ts = constant) for Tm − Ts = 1300◦C; it shows that the depth to the isotherms as defined by (2) are hyperbola. From this, one can readily see that if the thermal lithosphere is bounded by √ isotherms, the thickness of the lithosphere√increases as t. For back-of-theenvelope calculations you can use D ∼ 2.3 κt for lithospheric thickness. (For κ = 1 mm2 /s and t = 62.8 Ma, which is the average age of all ocean oceanic t (Myr) 0 0 50 100 150 200 y (km) 400 600 50 800 1000 100 150 Figure by MIT OCW. Figure 5.9: Oceanic geotherms. lithosphere currently at the Earth’s surface, D ∼ 104 km). This thickening occurs because the cool lithosphere reduces the temperature of the underlying material which can then become part of the plate. On the diagram the open circles depict estimates of lithospheric thickness from surface wave dispersion data. Note that even though the plate is moving and the resultant geometry is two dimensional the half space cooling model works for an observer that is moving along with the plate. Beneath this moving reference point the plate is getting thicker and thicker. Intermezzo 5.1 lithospheric thickness from surface wave dispersion In the seismology classes we have discussed how dispersion curves can be used to extract information about lithospheric structure from the seismic data. The thickness of the high wave speed ’lid’, the structure above the mantle low velocity zone, as determined from surface wave dispersion across parts of oceanic lithosphere of different age appears to plot roughly between the 900◦ and 1100◦ C isotherms, or at about T = 0.8 − 0.9Tm (see Figure 5.10). So the seismic ”lithosphere” seems to correspond roughly to thermal lithosphere. In other words; in short time scales — i.e. the time scale appropriate for seismic wave transmission (sec - min) — most of the thermal lithosphere may act as an elastic medium, whereas on the longer time scale the stress can be relaxed by steady state creep, in particular in the bottom half of the plate. However, a word of caution is in order since this interpretation of the dispersion data has been disputed. Anderson and co-workers argue (see, for instance, Anderson & Regan, GRL, vol. 10, pp. 183-186, 1983) that interpretation of surface wave dispersion assuming isotropic media results in a significant overestimation of the lid thickness. They have investigated the effects of seismic anisotropy and claim that the fast isotropic LID extends to a much cooler isotherm, at T ∼ 450◦ − 600◦ C, than the base of the thermal lithosphere. Heat flow If we know the temperature at the surface we can deduce the heat flow by calculating the temperature gradient: q = = ∂T ∂T ∂η = −k ∂z ∂η ∂z k(Tm − Ts ) ∂ k ∂ √ erf(η) − √ (Tm − Ts ) erf(η) = − ∂η ∂η 2 κt 2 κt −k (5.14) with ∂ erf(η) ∂η = ∂ 2 √ ∂η π �η 2 ∂ √ π ∂η �η 2 e−u du 0 = 2 2 2 e−u du = √ e−η π (5.15) 0 so that q = −k (Tm − Ts ) −η2 √ e πκt (5.16) 0 Age (Ma) 0 100 Depth (km) 300oC 50 900oC 600oC 1200oC 100 1325oC Figure by MIT OCW. Figure 5.10: Elastic thickness. For the heat flow proper we take z = 0 (q is measured at the surface!) so that η = 0 and q = −k 1 (Tm − Ts ) √ ⇒q∼ √ πκt t (5.17) with k the conductivity (do not confuse with κ, the diffusivity!). The important result is that according to the half-space cooling model the heat flow drops of as 1 over the square root of the age of the lithosphere. The heat flow can be measured, the lithospheric age t determined from, for instance, magnetic anomalies, and if we assume values for the conductivity and diffusivity, Eq. (5.17) can be used to determine the temperature difference between the top and the bottom of the plate √ πκt Tm − Ts = q (5.18) k Parsons & Sclater did this (JGR, 1977); assuming k = 3.13 JK−1 m−1 s−1 , CP = 1.17×103 Jkg−1 K−1 , and ρ = 3.3×103 kgm−3 and using q = 471 mWm−2 as the best fit to the data they found: Tm − Ts = 1350◦ ± 275◦ C. Comparison to observed heat flow data: qs (hfu) 6 4 2 0 0 50 100 150 t (106 yr) Figure by MIT OCW. Figure 5.11: . Near the ridge crest the observed heat flow is significantly lower than the heat flow predicted from the cooling half-space model. In old oceanic basins the heat flow seems to level off at around 46 mW/m2 , which suggests that beyond a certain age of the lithosphere the rate of conductive cooling either becomes smaller or the cooling is partly off set by additional heat production. Possible sources of heat which could prevent the half-space cooling are: 1. radioactivity (A is not zero!) 2. shear heating 3. small-scale convection below plate 4. hot upwelings (plumes) Intermezzo 5.2 Plate-cooling Models There are two basic models for the description of the cooling oceanic lithosphere, a cooling of a uniform half space and the cooling of a layer with some finite thickness. The former is referred to as the half-space model (first described in this context by Turcotte and Oxburgh, 1967); The latter is also known as the plate model (first described by McKenzie, 1967). Both models assume that the plate moves as a unit, that the surface of the lithosphere is at an isothermal condition of 0◦ C, and that the main method of heat transfer is conduction (a good assumption, except at the ridge crest). The major difference (apart from the mathematical description) is that in the halfspace model the base of the lithosphere is defined by an isotherm (for instance 1300◦ C) so that plate thickness can grow indefinitely whereas in the plate model the plate thickness is limited by some thickness L. The two models give the same results for young plates near the ridge crest, i.e. the thickness is such that the ”bottom” of the lithosphere is not yet ”sensed”. However, they differ significantly after 50 Myr for heat flow predictions and 70 Myr for topography predictions. It was realized early on that at large distances from the ridge (i.e., large ages of the lithosphere) the oceans were not as deep and heat flow not as low as expected from the half-space cooling model (there does not seem to be much thermal difference between lithosphere of 80 and 160 Myr of age). The plate model was proposed to get a better fit to the data, but its conceptual disadvantage is that it does not explain why the lithosphere has a maximum thickness of L. The half space model makes more sense physically and its mathematical description is more straightforward. Therefore, we will discuss only the half space cooling model, but we will also give some relevant comparisons with the plate model. 5.4 Thermal structure of the oceanic lithosphere Bathymetry The second thermal effect on the evolution of the cooling lithosphere is its subsidence or the increase in ocean depth with increasing age. This happens because when the mantle material cools and solidifies after melting at the MOR it is heavier than the density of the underlying mantle. Since we have seen that the plate thickens with increasing distance from the MOR and if the plate is not allowed to subside this would result in the increase in hydrostatic pressure at some reference depth. In other words the plate would not be in hydrostatic equilibrium. But when the lithosphere subsides, denser material will be replaced by lighter water so that the total weight of a certain column remains the same. The requirement of hydrostatic equilibrium gives us the lateral variation in depth to the ocean floor. Application of the isostatic principle gives us the correct ocean floor topography. Let ρw and ρm represent the density of water (ρw ∼ 1000 kg/m3 ) and the mantle/asthenosphere (ρm ∼ 3300 kg/m3 ), respectively, and ρ(z, t) the density as a function of time and depth within the cooling plate. The system is in hydrostatic equilibrium when the total hydrostatic pressure of a column under the ridge crest at depth w + zL is the same as the pressure of a column of the x2 x1 x Ridge w �w WATER � yL LITHOSPHERE �m |x | t2 = u2 |x | t1 = u1 y ASTHENOSPHERE Figure by MIT OCW. Figure 5.12: Oceanic isostasy. same width at any distance from the MOR: �zL ρ(z, t) dz (5.19) �zL w(ρw − ρm ) + (ρ − ρm ) dz = 0 (5.20) ρm (w + zL ) = wρw + 0 or 0 According to (5.20) the mass deficiency caused by w(ρw − ρm ) (which is less than 0!) is balanced by the difference between ρ − ρm (>0!) integrated over the (as yet unknown) lithospheric thickness. The lithosphere is thus heavier than the underlying half space! (Assuming, as we do here, that the lithosphere has the same composition as the asthenosphere). This increase in density is due to cooling; the relationship between the change in density due to a change in temperature is given by dρ = −αρdT ⇒ ρ − ρm = −αρ(T − Tm ) ⇒ ρ(z, t) = ρm + αρ(Tm − T (z, t)) (5.21) (5.22) with α the coefficient of thermal expansion, so that �zL w(ρw − ρm ) + ρm α (Tm − T ) dz = 0 0 With T = T (z, t), this gives (verify!!) (5.23) � �� �zL � z dz w(ρw − ρm ) = ρm α (Tm − Ts ) − (Tm − Ts ) erf √ 2 κt 0 � �� �zL � z = ρm α(Tm − Ts ) 1 − erf √ dz 2 κt 0 = ρm α(Tm − Ts ) �zL erfc � 0 z √ 2 κt � dz (5.24) we can change the integration boundary from zL to ∞ because at the base of the lithosphere T → Tm and ρ → ρm so that we can take the compensation a t any depth beneath the base of the cooling lithosphere (and erfc integrated √ from zL to ∞ is very small). If we also use η = z/(2 κπ) (see above) than we can write w(ρw − ρm ) = ρm α(Tm − Ts ) �∞ 0 � z erfc √ 2 κt � dz ∞ √ � = 2ρm α(Tm − Ts ) κt erfc(η) (5.25) 0 now use �∞ 0 w(t) = erfc(q) dq = √1 π 2ρm α(Tm − Ts ) (ρw − ρm ) to get � κt π (5.26) with w(t) the depth below the ridge crest; if the crest is at depth w0 (5.26) becomes � 2ρm α(Tm − Ts ) κt w(t) = w0 + (5.27) (ρw − ρm ) π So from the half-space cooling model it follows that the depth to the sea floor Sclater increases as the square root of age! Using α = 3.2×10−5◦ C−1 , Pars ons & √ (1977) found from the fit to bathymetry data gives: w(t) = 2500 + +350 t [m], for t < 70 Ma. There is still a lively debate about the details of the parameters that give the best fit to the model, see, for instance the papers by Stein & Stein (Nature, 1992) and McNutt (Reviews of Geophysics, 1995). But despite the ongoing discussions it is fair to say that these models have been very successful in predicting heat flow, topography, gravity, and have thus played a major role in the understanding of the evolution of oceanic lithosphere with time. The typical game that is played by such successful theoretical models, in particular ones that are so simple (= easy + fast to compute) as the cooling models is to predict the first order behavior of a certain process and take out that trend from the observed data. In this case, the residual signal is then analyzed for deviations from the simple conduction model. The addition of heat to the system (for instance by plumes) could cause anomalous topography (”thermal topography”) whereas the effect of deep dynamic processes in the Earth’s mantle can cause ”dynamic topography”. Removing the effects of conduction alone thus helps to isolate the structural signal due to other processes. This is likely to continue, perhaps with the new model (GHD1) by Stein & Stein (1992) instead of that by earlier workers; since regional differences are often larger than the residual between observed and predicted heat flow and depth curves one could question how useful a (set of) simple model(s) is (are). For instance, if one allows the thermal expansion coefficient as a free parameters in the inversions, one might also look into allowing lateral variation of this coefficient. Davies and Richards argue that the success of the cooling models in predicting the topography and heat flow over almost the entire age range of oceanic lithosphere (they attribute the deviations to the choice of the wrong sites for data — which is rather questionable) indicates that conductive cooling is the predominant mode of heat loss of most of the lithosphere (about 85% of the heat lost from the mantle flows through oceanic lithosphere), which suggests that the lithosphere is the boundary layer of a convective system with a typical scale length defined by the plates (plate-scale flow). They follow up on a concept tossed up by Brad Hager that the oceanic lithosphere organizes the flow in the deeper mantle. It is for arguments such as these that question as to whether or not the topography levels off after, say, 80Ma, in not merely of interest to statisticians. It is quite clear that the details of the bathymetry of the oceans still contain significant keys to the understanding of dynamic processes in the deep interior of the Earth. One of the remarkable aspects of the square root of time variation of ocean depth is that it does a very good job in describing the true bathymetry, even at distances pretty close to the MORs. This indicates that conduction alone is likely to be the predominant mode of heat loss, even close to the MOR. The absence of any substantial dynamic topography near the ridge crest suggests that the active, convection related upwellings are not significant. The upwelling is ”passive”: the plates are pulled apart (mainly as the result of the gravitational force, the slab pull, acting on the subducting slabs) and the asthenospheric mantle beneath the ridges flows to shallower depth to fill the vacancy. In doing so the material will cross the solidus, the temperature at which rock melts (which decreases with decreasing depth) so that the material melts. This process is known as decompression melting (see Turcotte & Schubert, Chapter 1), which results in a shallow magma chamber beneath the MOR instead of a very deep plume-like conduit. Temperature (oC) Depth (km) 1000 0 50 1200 1400 Depth at which partial melting begins Temperature of ascending mantle rock Solidus 100 Figure by MIT OCW. Figure 5.13: Pressure-release melting. 5.5 Bending, or flexure, of thin elastic plate Introduction We have seen that upon rifting away from the MOR the lithosphere thickens (the base of the thermal lithosphere is defined by an isotherm, usually Tm ≈ 1300◦C) and subsides, and that the cooled lithosphere is more dense than the underlying mantle. In other words, it forms a gravitationally unstable layer. Why does it stay atop the asthenosphere instead of sinking down to produce a more stable density stratification? That is because upon cooling the lithosphere also acquires strength. Its weight is supported by its strength; the lithosphere can sustain large stresses before it breaks. The initiation of subduction is therefore less trivial than one might think and our understanding of this process is still far from complete. The strength of the lithosphere has important implications: 1. it means that the lithosphere can support loads, for instance by seamounts 2. the lithosphere, at least the top half of it, is seismogenic 3. lithosphere does not simply sink into the mantle at trenches, but it bends or flexes, so that it influences the style of deformation along convergent plate boundaries. Investigation of the bending or flexure of the plate provides important information about the mechanical properties of the lithospheric plate. We will see that the nature of the bending is largely dependent on the flexural rigidity, which in turn depends on the elastic parameters of the lithosphere and on the elastic thickness of the plate. An important aspect of the derivations given below is that the thickness of the elastic lithosphere can often be determined from surprisingly simple observations and without knowledge of the actual load. In addition, we will see that if the bending of the lithosphere is relatively small the entire mechanical lithosphere behaves as an elastic plate; if the bending is large some of the deformation takes place by means of ductile creep and the part of the lithosphere that behaves elastically is thinner than the mechanical lithosphere proper. 5.5.1 Basic theory To derive the equations for the bending of a thin elastic plate we need to 1. apply laws for equilibrium: sum�of the forces is zero and the sum of all � moments is zero: F = 0 and M =0 2. define the constitutive relations between applied stress σ and resultant strain ǫ 3. assume that the deflection w ≪ L, the typical length scale of the system, and h, the thickness of the elastic plate ≪ L. The latter criterion (#3) is to justify the use of linear elasticity. Figure 5.14: Deflection of a plate under a load. In a 2D situation, i.e., there is no change in the direction of y, the bending of a homogeneous, elastic plate due to a load V (x) can be described by the fourth-order differential equation that is well known in elastic beam theory in engineering: D d4 w d2 w + P 2 = V (x) 4 dx dx (5.28) with w = w(x) the deflection, i.e., the vertical displacement of the plate, which is, in fact, the ocean depth(!), D the flexural rigidity, and P a horizontal force. The flexural rigidity depends on elastic parameters of the plate as well as on the thickness of the plate: D= Eh3 12(1 − ν 2 ) (5.29) with E the Young’s modulus and ν the Poisson’s ratio, which depend on the elastic moduli µ and λ (See Fowler, Appendix 2). The bending of the plate results in bending (or fiber) stresses within the plate, σxx ; depending on how the plate is bent, one half of the plate will be in compression while the other half is in extension. In the center of the plate the stress goes to zero; this defines the neutral line or plane. If the bending is not too large, the stress will increase linearly with increasing distance z ′ away from the neutral line and reaches a maximum at z ′ = ±h/2. The bending stress is also dependent on the elastic properties of the plate and on how much the plate is bent; σxx ∼ elastic moduli ×z ′ × curvature, with the curvature defined as the (negative of the) change in the slope d/dx(dw/dx): σxx = − Eh3 ′ d2 w z 1 − ν 2 dx2 (5.30) y αxx M comp h y=0 x M ext αxx Figure by MIT OCW. Figure 5.15: Curvature of an elastic plate. This stress is important to understand where the plate may break (seismicity!) with normal faulting above and reverse faulting beneath the neutral line. The integrated effect of the bending stress is the bending moment M , which results in the rotation of the plate, or a plate segment, in the x − z plane. h M= �− 2 σxx z ′ dz ′ (5.31) h 2 Equation (5.28) is generally applicable to problems involving the bending of a thin elastic plate. It plays a fundamental role in the study of such problems as the folding of geologic strata, the development of sedimentary basins, the post-glacial rebound, the proper modeling of isostasy, and in the understanding of seismicity. In class we will look at two important cases: (1) loading by sea mounts, and (2) bending at the trench. Before we can do this we have to look a bit more carefully at the dynamics of the system. If we apply bending theory to study lithospheric flexure we have to realize that if some load V or moment M causes a deflection of the plate there will be a hydrostatic restoring force owing to the replacement of heavy mantle material by lighter water or crustal rock. The magnitude of the restoring force can easily be found by applying the isostasy principle and the effective load is thus the applied load minus the restoring force (all per unit length in the y direction): V = Vapplied − Δρwg with w the deflection and g the gravitational acceleration. This formulation also makes clear that lithospheric flexure is in fact a compensation mechanism for isostasy! For oceanic lithosphere Δρ = ρm − ρw and for continental flexure Δρ = ρm − ρc . The bending equation that we will consider is thus: D d4 w d2 w + P 2 + Δρwg = V (x) 4 dx dx (5.32) Loading by sea mounts Let’s assume a line load in the form of a chain of sea mounts, for example Hawaii. Figure 5.16: Deflection of an elastic plate under a line load. Let V0 be the load applied at x = 0 and V (x) = 0 for x �= 0. With this approximation we can solve the homogeneous form of (5.32) for x > 0 and take the mirror image to get the deflection w(x) for x < 0. If we also ignore the horizontal applied force P we have to solve D d4 w + Δρwg = 0 dx4 (5.33) The general solution of (5.33) is x w(x) = e α � � x x x� x x� A cos + B sin + e− α C cos + D sin α α α α (5.34) with α the flexural parameter, which plays a central role in the extraction of structural information from the observed data: α= � 4D Δρg � 14 (5.35) Courtesy of Annual Reviews Inc. Used with permission. Figure 5.17: . The constants A − D can be determined from the boundary conditions. In this case we can apply the general requirement that w(x) → 0 for x → ∞ so that A = B = 0, and we also require that the plate be horizontal directly beneath x = 0: dw/dx = 0 for x = 0 so that C = D: the solution becomes x w(x) = Ce− α � cos x x� + sin α α (5.36) From this we can now begin to see the power of this method. The deflection w as a function of distance is an oscillation with period x/α and with an exponentially decaying amplitude. This indicates that we can determine α directly from observed bathymetry profiles w(x), and from equations (5.36) and (5.29) we can determine the elastic thickness h under the assumption of values for the elastic parameters (Young’s modulus and Poisson’s ratio). The flexural parameter α has a dimension of distance, and defines, in fact, a typical length scale of the deflection (as a function of the “strength” of the plate). The constant C can be determined from the deflection at x = 0 and it can be shown (Turcotte & Schubert) that C = (V0 α3 )/(8D) ≡ w0 , the deflection x/� 0 1 2 3 w/w0 x0/� 4 5 xb/� 0.5 1 Figure by MIT OCW. Figure 5.18: A deflection profile. beneath the center of the load. The final expression for the deflection due to a line load is then V0 α3 − x � x x� x≥0 (5.37) e α cos + sin w(x) = 8D α α Let’s now look at a few properties of the solution: • The half-width of the depression can be found by solving for w = 0. From (5.37) it follows that cos(x0 /α) = − sin(x0 /α) or x0 /α = tan−1 (−1) ⇒ x0 = α(3π/4 + nπ), n = 0, 1, 2, 3 . . . For n = 0 the half-width of the depression is found to be α3π/4. • The height, wb , and location, xb , of forebulge ⇒ find the optima of the solution (5.37). By solving dw/dx = 0 we find that sin(x/α) must be zero ⇒ x = nπα, and for those optima w = w0 e−nπ , n = 0, 1, 2, 3 . . . For the location of the forebulge: n = 1, xb = πα and the height of the forebulge wb = −w0 e−π or wb = −0.04w0 (very small!). Important implications: The flexural parameter can be determined from the location of either the zero crossing or the location of the forebulge. No need to know the magnitude of the load! The depression is narrow for small α, which means either a weak plate or a small elastic thickness (or both); for a plate with large elastic thickness, or with a large rigidity the depression is very wide. In the limit of very large D the depression is infinitely wide but the amplitude w0 , is zero ⇒ no depression at all! Once α is known, information about the central load can be obtained from Eq. (5.37) Note: the actual situation can be complicated by lateral variations in thickness h, fracturing of the lithosphere (which influences D), compositional layering within the elastic lithosphere, and by the fact that loads have a finite dimensions. Flexure at deep sea trench With increasing distance from the MOR, or with increasing time since formation at the MOR, the oceanic lithosphere becomes increasingly more dense and if the conditions are right5 this gravitational instability results in the subduction of the old oceanic plate. The gravitational instability is significant for lithospheric ages of about 70 Ma and more. We will consider here the situation after subduction itself has been established; in general the plate will not just sink vertically into the mantlebut it will bend into the trench region. -V0 wb x=0 x -M0 x = x0 x = xb Figure by MIT OCW. Figure 5.19: This bending is largely due to the gravitational force due to the negative buoyancy of the part of the slab that is already subducted M0 . For our modeling 5 Even for old oceanic lithosphere the stresses caused by the increasing negative buoyancy of the plate are not large enough to break the plate and initiate subduction. The actual cause of subduction initiation is still not well understood, but the presence of pre-existing zones of weakness (e.g. a fracture zone, thinned lithosphere due to magmatic activity — e.g. an island arc) or the initiation of bending by means of sediment loading have all been proposed (and investigated) as explanation for the triggering of subduction. we assume that the bending is due to an end load V0 and a bending moment M0 applied at the tip of the plate. As a result of the bending moment the slope � 0 at x = 0 (note the difference with the seamount example where this dw/dx = slope was set to zero!). The important outcome is, again, that the parameter of our interest, the elastic thickness h, can be determined from the shape of the plate, in vertical cross section, i.e. from the bathymetry profile w(x)!, in the subduction zone region, without having to know the magnitudes of V0 and M0 . We can use the same basic equation (5.33) and the general solution (5.34) (with A = B = 0 for the reason given above) � x x x� w(x) = e− α C cos + D sin (5.38) α α but the boundary conditions are different and so are the constants C and D. At x = 0 the bending moment6 is −M0 and the end load −V0 . It can be shown (Turcotte & Schubert) that the expressions for C and D are given by C = (V0 α + M0 ) α2 2D and D = − M0 α2 2D (5.39) so that the solution for bending due to an end load and an applied bending moment can be written as α2 e−x/α � x x� w(x) = (V0 α + M0 ) cos − M0 sin (5.40) 2D α α We proceed as above to find the locations of the first zero crossing and the fore bulge, or outer rise. w(x) = 0 ⇒ tan(x0 /α) = 1 + αV0 /M0 dw/dx = 0 ⇒ tan(xb /α) = −1 − 2M0 /αV0 (5.41) (5.42) In contrast to similar solutions for the sea mount loading case, these expressions for x0 and xb still depend on V0 and M0 . In general V0 and M0 are unknown. They can, however, be eliminated, and we can show the dependence of w(x) on x0 and xb , which can both be estimated from the nathymetry profile. A perhaps less obvious but elegant way of doing this is to work out tan(1/α(x0 − xb )). Using sine and cosine rules (see Turcotte & Schubert, 3.17) one finds that � � xb − x0 tan =1 (5.43) α so that x0 − xb = (π/4 + nπ)α, n = 0, 1, 2, 3, . . . For n = 0 one finds that α = 4(x0 − xb )/π, so that the elastic thickness h can be determined if one can measure the horizontal distance between x0 and xb . 6 At this moment, it is important that you go back to the original derivation of the plate equation in Turcotte & Schubert and realize they obtained their results with definite choices as to the signs of applied loads and moments — hence the negative signs. 2 1 100 -100 x (km) 200 300 w (km) -1 -2 -3 -4 Figure by MIT OCW. Figure 5.20: After a bit of algebra one can also eliminate α to find the deflection w(x) as a function of wb , x0 , and xb . The normalized deflection w/wb as a function of normalized distance (x − x0 )/(xb − x0 ) is known as the Universal Flexure Profile. � � �� � � �� π x − x0 π x − x0 w(x) √ π sin − (5.44) = 2e 4 exp − 4 xb − x0 wb 4 xb − x0 In other words, there is a unique way to bend a laterally homogeneous elastic plate so that it goes through the two points (x0 , 0) and (xb , wb ) with the condition that the slope is zero at x = xb . The example of the Mariana trench shown in Figure 5.20 demonstrates the excellent fit between the observed bathymetry and the prediction after Eq. (5.44) (for a best fitting elastic thickness h as determined from the flexural parameter calculated from equation (5.43). Bending stress and seismicity Many shallow earthquakes occur in near the convergent margin. Both in the overriding plates as well as in the subducting plate. The latter can be attributed to the bending stresses in the plate. The bending stress is given by Eq. (5.30). Earthquakes are most likely to occur in the region where the bending stress is largest (that’s the place where the elastic plate is most likely to fail if there are no pre-existing inhomogeneities such as transform faults). To find the horizontal location where the stress is largest we must solve � � dσxx d d2 w d3 w =0 ⇒ (5.45) = 0 = dx dx dx2 dx3 This gives the location x where the stress is a maximum (or minimum!) and substitution in (5.30), with the flexural parameter D determined as above from the bathymetry profile, then gives the amplitude of the maximum stress. If this stress exceeds the strength of the plate, failure will occur. The mechanism of the earthquake depends on the location relative to the neutral stress plane. 5.6 The upper mantle transition zone Derivation of density variation with depth : Adams-Williamson Equation How about density? Can the radial variation in density and the elastic moduli be constrained independently from the travel time curves? Indirectly, yes! This was first shown by Adams and Williamson in 1923. Here, we will only give the basic principles and, in particular, discuss its implications for our understanding of the Earth’s physical state. The fundamental result I want you to remember is that the Adams-Williamson equation relates the gradient in density to radial variations in seismic wave speed (through the seismic parameter) and the mass of the Earth, which quantities are assumed to be known, but that this result only applies to homogeneous regions of the same physical phase. From the travel time curves we can determine radial variations of P and S-wave speed, α(r) and β(r). α2 = β2 = κ + 4/3µ ρ µ ρ (5.46) (5.47) which can be combined to get what is known as the seismic parameter 4 κ Φ = α2 − β 2 = 3 ρ (5.48) where α, β, ρ, µ, κ, and Φ are all functions of radius:α(r), β(r), ρ(r), µ(r), κ(r), and Φ(r). The seismic parameter is also known as the bulk sound velocity, as the counterpart of the shear velocity β. (Notice that the incompressibility κ in these equations is, in fact, the adiabatic incompressibility or bulk modulus κS because the time scale of any change in κ due to changes in temperature T are much larger than the transit time of a seismic wave.) The aim is to show that not only the density-normalized shear and bulk moduli can be determined, but also the density itself (and thus µ and κ). In general, variations in density can be due to changes in pressure (dP ), temperature (dT ), composition (dc) and physical phase (dϕ), which can be written (in gradient form) as: dρ = dr � ∂ρ ∂P � dP + dr � ∂ρ ∂T � dT + dr � ∂ρ ∂c � dc + dr � ∂ρ ∂ϕ � dϕ dr (5.49) For a homogeneous medium (same composition and phase throughout) this equation simplifies to; dρ = dr � ∂ρ ∂P � dP + dr � ∂ρ ∂T � dT dr (5.50) For the sake of the argument I will concentrate on the effect of adiabatic compression, i.e., there is no variation of density with temperature. dρ ≈ dr � ∂ρ ∂P � dP dr (5.51) This assumption seems reasonable for most of the convecting mantle, and leads to the original Adams-Williamson equation. For thermal boundary layers such as the lithosphere and the lowermost mantle (D”), and —- in case of layered convection — a TBL between the upper and the lower mantle, an additional gradient term has to be taken into account, and this modification has been applied by Birch in his famous 1952 paper (see Fowler §4.3, and Stacey §5.3.1). For adiabatic self-compression the increase in pressure that results from the descent from radius r + dr to radius r is due to the weight of the overlying shell with thickness δr, so that the pressure gradient can be written as: dP = −gρ, dr with g=G m(r) r2 (5.52) The other term in Eq. (5.51), the pressure derivative of the density, can be evaluated in terms of the adiabatic bulk modulus κS : κS = dP dP increase in pressure =− =ρ dV /V fractional change in volume dρ (5.53) Substitution of (5.52) and (5.53) in (5.51) and using (5.48) gives the AdamsWilliamson equation: dρ =− dr � ρ κS � ρGm(r) ρGm(r) =− r2 Φr2 (5.54) which relates the density gradient to the known seismic parameter and the gravitational attraction of the mass m(r). Rewrite for m(r) m(r) = 4π � 0 r 2 ρ(a)a da = MEarth − 4π � REarth ρ(a)a2 da (5.55) r shows that m(r) is, in fact, the mass of the Earth less the mass of the shell between point r and the radius of the Earth REarth . The mass of the Earth is assumed to be know from astronomical data and is an important constraint on the density gradient. So the only unknown in (5.55) is the density ρ(a) between r and REarth . We can find a solution of (5.54) by working from the Earth’s surface to larger depths: at the surface, the density of crustal rock is fairly well known so that the density gradient can be determined for the crust. This gradient is then used to estimate the density at the base of the crust, which is then used to calculate the mass of the crustal shell. In this way we can carry on the differentiation and integration to larger depths. As already mentioned, and explicit in (5.55), any solution of (5.54) must agree with the total mass of the Earth, as well as with the moment of inertia, which forms the second independent constraint on the solution. This process can only be applied in regions where dα/dr = dϕ/dr = 0, and in this form one must also require ∂ρ/∂T = 0, but — as mentioned above — there are approximations to (5.54) that take small deviations from adiabatic compression into account. Application of the (modified) Adams-Williamson equation by, amongst others, Bullen resulted in pretty good density models for the Earth. The upper mantle transition zone In 1952, Birch realized that both the density gradient and the wave speed gradients in the Earth mantle between 200 and 900 km in depth are larger than expected from adiabatic compression only, see the abstract of his famous paper (attached). This means that either dc/dr �= 0 or dφ/dr �= 0, or both. The mantle region where the density and wave speed gradient are larger than predicted from adiabatic compression alone is loosely referred to as the (upper mantle) transition zone7 . There is no consensus yet on which one applies to the Earth but it is now clear that dϕ/dr = � 0 is a sufficient condition and is probably the most important factor to explain the excess density. Birch suggested that phase transformations in the Mg, Fe silicate system (Mg, Fe)2 SiO4 (olivine, spinel) and (Mg, Fe)SiO3 (pyroxene) could explain the increase in density required to explain the nonadiabatic parts of the density and wave speed gradients. We now know that phase transformations do indeed occur in this mineral system at depths of about 410 and 660 km. Initially the sharpness of the interface as deduced from the 7 This is the original definition of the transition zone. Later, it became common to use the term “transition zone” in a more restricted to indicate the mantle region between the 410 and 660 km discontinuities. In terms of mantle processes (convection, slab behavior) the original definition is more useful. Density (gcm-3) 0 3 4 5 6 7 8 9 10 11 12 13 1000 Depth (km) 2000 3000 4000 5000 6000 Figure by MIT OCW. Figure 5.21: The density of the Earth according to model ak135. reflection and phase conversion of high frequency seismic waves was used as evidence for compositional layering and against effects of a phase change. However, experimental rock mechanics in the late eighties demonstrated that phase changes can occur over sufficiently narrow depth ranges to explain the seismic observations, see attached phase diagrams by Ito and Takahashi (JGR, 1989). The phase changes in the (Mg, Fe) silicates play an important role in mantle dynamics because the pressure induced phase changes are also temperature dependent! This means that phase changes can occur at different depth depending on the temperature of the medium in which the transformation occurs. The temperature dependence is governed by the value of the Clapeyron slope dP/dT of the boundaries between the stability fields of Olivine (Ol), Spinel (Sp), and Perovskite/Magnesiowüstite (Pv+Mw) in the P − T diagram. The phase diagrams by Ito & Takahashi at 1100◦C and 1600◦C illustrate that the phase Figure 5.22: Phase diagrams of mineral transformations in the mantle. change occurs at smaller pressure if the temperature increase; i.e., the Clapeyron slope for the transition from Sp → Pv+Mw is negative! It’s a so called endothermic phase change: upon phase transformation the material looses heat and cools down. In contrast, Ol → Sp transition that marks the phase change at about 410 km depth has a positive Clapeyron slope and is exothermic, i.e. there is a release of latent heat upon transformation and the material will warm up. What does this mean for dynamics and plate driving forces? In the diagram I have given schematically the stability fields of Ol, Sp, and Pv+Mw, and the boundaries between them (i.e. the Clapeyron slopes). If one would descend into the mantle along an average mantle geotherm one would cross the Clapeyron slope where Ol and Sp coexist at a pressure that corresponds to a depth of about 410 km and the phase line between Sp and Pv+Mw at a pressure corresponding to about 660 km depth. Consider now the situation that a slab of cold, former oceanic lithosphere subducts into the mantle and crosses the stability fields of the silicates (assume for simplicity that the composition of the slab is the same as the mantle — which is not the case). Within the slab the phase transformation from Ol → Sp will occur at a shallower depth than in the ambient mantle. This means that for depths just less than 410 km the density within the slab is locally higher than in the ambient mantle, and this, in fact, gives rise to an extra negative buoyancy force that helps the slab to subduct (it is an important plate driving force). At 660 km the dynamical effect of the phase change is the opposite. Inside and in the direct vicinity of the slab the phase boundary will be depressed; consequently, the density in the slab is less than the density of the ambient mantle which creates a buoyancy force that will resist the further penetration of the slab. Figure 5.23: Effects of phase transformations on downgoing slabs. From the diagram it is clear that the steeper the Clapeyron slope, the stronger the dynamic effects. On the one hand, a lot of laboratory research is focused on estimating the slopes of these phase boundaries in experimental conditions. On the other hand, and that brings us back to seismology, seismologists attempt to estimate the topography of the seismic discontinuity and thus constrain the clapeyron slope and asses the dynamical implications. Important classes of seismological data that have the potential to constrain both the sharpness of and the depth to the discontinuities are reflections and phase (mode) transformations. An example of a useful reflection is the underside reflection of the P KIKP P KIKP , or P KPDF P KPDF , or just P ′ P ’, at the 660 km discontinuity. ′ Since the paths of the P ′ P ′ phase and the P660 P ′ are almost similar except for near the reflection point, the difference in travel time gives direct information about the depth to the interface. Another example is the use of SS underside reflections S660 S. Apart from proper phase identification (usually one applies stacking techniques to suppress signal to noise) the major problem with such techniques is that one has to make assumptions about upper mantle structure between the Earth’s surface and the discontinuity, and these corrections are not always reliable. The time difference between the reflections at the surface and the discontinuity contains information about the depth to the interface, whereas the frequency content of both the direct and the reflected phase gives information about the sharpness of the interface. Also phase conversions can be used! This line of research is still very active, and there is some consensus only about the very long wave length variations in depth to the seismic discontinuities.