Dynamic origins of seismic wavespeed variation in D′′

The D′′ discontinuity is defined by a seismic velocity increase of 1–3% about 250 km above the core-mantle boundary (CMB), and is mainly detected beneath locations of inferred paleosubduction. A phase change origin for the interface can explain a triplicated arrival observed in seismic waveform data and is supported by the recent discovery of a post-perovskite phase transition. We investigate the interaction of slabs, plumes, and the phase change within D′′ in 2-D compressible convection calculations, and predict waveform complexity in synthetic seismic data. The dynamic models produce significant thermal and phase heterogeneity in D′′ over small distances and reveal a variety of behaviors including: (1) distinct pPv blocks separated by upwellings, (2) notches at the top of a pPv layer caused by plume heads, (3) regions of Pv embedded within a pPv layer due to upwellings. Advected isotherms produce complicated thermal structure that enable multiple crossings of the phase boundary. Perturbations to S, SdS, and ScS arrivals (distances < 84 degrees) are linked to the evolutionary stage of slabs and plumes, and can be used to determine phase boundary height and velocity increase, volumetric wavespeed anomaly beneath the discontinuity, and possibly the lengthscale of slab folding near the CMB. Resolving fine-scale structure beneath the interface requires additional seismic phases (e.g., Sd, SKS) and larger distances (> 80 degrees).


Introduction
The D discontinuity is characterized by a seismic velocity increase of 1-3% approximately 250 kilometers above the core-mantle boundary (CMB) (see review by Wysession et al., 1998). Seismic waveform modeling detects a velocity jump beneath locations of inferred paleosubduction, including Alaska, the Caribbean, Central America, India, and Siberia. Furthermore, in seismic tomography the high velocity and deep anomalies in these regions are interpreted as slabs (e.g., Grand, 2002).
The discontinuity can explain a triplicated arrival SdS (PdP) between S (P) and ScS (PcP) observed in waveform data (Lay and Helmberger, 1983). Between epicentral distances 65-83 degrees, SdS is a composite arrival from a discontinuity reflection Sbc (Pbc) and a ray turning below the interface Scd (Pcd). The horizontally-polarized S-wave (SH) triplicated arrival is often analyzed at ∼ 10 s period. SH is usually not contaminated by other phases or mode conversions and modern array processing techniques using broadband data can recognise energy on the tangential component that originates from other phenomena such as SKS splitting (e.g., Garnero and Lay, 2003;Wookey and Kendall, 2007). Shorter periods (∼ 1 s) can detect the P-wave triplication (PdP) but data stacking is required to suppress noise.
There are generally fewer or weaker SdS (PdP) detections outside highvelocity areas which suggests subduction history influences the presence and strength of the arrival. In contrast to seismic observations, a pre-existing basal chemical layer generates strong (weak) SdS below upwellings (downwellings)  and does not generate short wavelength heterogeneity on the discontinuity (Tackley, 1998). However, detections of SdS beneath the seismically slow central Pacific suggest the discontinuity height above the CMB or velocity increase may be modulated by composition in some regions (Garnero et al., 1993;Lay et al., 2006).
Early dynamic calculations and waveform modeling propose a thermal slab interacting with a phase change with a positive Clapeyron slope can explain the origin of the D discontinuity (Sidorin et al., 1999a. Incident rays are refracted by a high velocity thermal slab above D and turn beneath the discontinuity in the higher velocity phase to produce the SdS arrival. Furthermore, a high-temperature thermal boundary layer generates a negative seismic velocity gradient at the CMB. This counteracts the velocity step to ensure that differential travel times (e.g., ScS-S) and diffracted waveforms match data Lay, 1987, 1990). Sidorin et al. (1999b) unify regional seismic models from waveform studies by proposing a global discontinuity with an ambient phase change elevation of 200 km above the CMB and a Clapeyron slope of 6 MPa K −1 . The actual boundary locally is thermally perturbed upward and downward.
The subsequent discovery of a phase transition in MgSiO 3 from silicate perovskite (Pv) to "post-perovskite" (pPv) at CMB conditions supports the phase change hypothesis Oganov and Ono, 2004). Ab initio calculations at 0 K suggest that the transformation of isotropic aggregates of Pv to pPv increases the S-wave velocity by 1 to 1.5% and changes the P-wave velocity by -0.1 to 0.3% (Oganov and Ono, 2004;Tsuchiya et al., 2004;Iitaka et al., 2004). First principle calculations at high temperature produce similar wavespeed variations Wentzcovitch et al., 2006), and high pressure experiments resolve a 0-0.5% increase in S-wave velocity (Murakami et al., 2007).
Lattice-preferred orientation in the pPv phase may reconcile the small variation in wavespeeds from mineral physics with the 1-3% increase across the D discontinuity deduced from seismic data (e.g., Wookey et al., 2005). Fe content in Pv may decrease (Mao et al., 2004) or increase (Tateno et al., 2007) the phase transition pressure, and Al may broaden the mixed phase region (Catalli et al., 2010;Akber-Knutson et al., 2005). However, other researchers find that compositional variations have little effect (Hirose et al., 2006;Murakami et al., 2005). Independent of the width of the mixed phase region, strain-partitioning into weak pPv can produce a sharp seismic discontinuity due to the rapid change in seismic anisotropy (Ammann et al., 2010). Thomas et al. (2011) reconcile a negative P-wave contrast from D beneath the Caribbean with a positive P-wave contrast beneath Eurasia by proposing a Pv-pPv phase transition with a fraction of 12% alignment in pPv.
The phase transformation destabilizes the lower thermal boundary layer, which produces more frequent upwellings (Sidorin et al., 1999a;Nakagawa and Tackley, 2004). A dense basal layer interacting with a chemically heterogeneous slab and the Pv-pPv transition produces strong lateral gradients in composition and temperature that may also explain finer structure within D (Tackley, 2011). A "double crossing" of the phase boundary occurs when pPv transforms back to Pv just above the CMB and may explain neighboring seismic discontinuities (Hernlund et al., 2005). However, the precritical reflection from the second putative phase transition is a relatively low amplitude arrival and therefore difficult to detect (Flores and Lay, 2005;Sun et al., 2006). The Clapeyron slope and transition temperature (or pressure) control the topography of the pPv layer (e.g., Monnereau and Yuen, 2007), The CMB beneath Central America is illuminated through the propagation of seismic waves along a narrow corridor from South American subduction zone events to broad-band networks in North America. Furthermore, seismic tomography and plate reconstructions indicate deeply penetrating slab material (e.g., Grand, 2002;Ren et al., 2007). Under the Cocos Plate, the D discontinuity can be modeled at constant height above the CMB with S-wave variations of 0.9-3.0% Ding and Helmberger, 1997) or as an undulating north-south dipping structure from 300 km to 150 km above the CMB with constant D velocity (Thomas et al., 2004). The Pv-pPv transition may account for the positive jump seen in S-wave models and smaller negative P-wave contrast (Hutko et al., 2008;Kito et al., 2007;Wookey et al., 2005). Furthermore, its interaction with a buckled slab may explain an abrupt 100 km step in the discontinuity (Sun and Helmberger, 2008;Hutko et al., 2006;Sun et al., 2006). A second deeper negative reflector appears in some locations about 200 km below the main discontinuity. This may originate from the base of a slab, a plume forming beneath a slab (Tan et al., 2002), back-transformation of pPv-Pv, chemical layering (Kito et al., 2007;Thomas et al., 2004), or out-of-plane scatterers (Hutko et al., 2006).
It is now timely to reanalyze the role of slabs in the deep mantle and their seismic signature following discovery of the Pv-pPv transformation and a proposed "double crossing" of the phase boundary. In this study we follow a similar approach to  using a compressible formulation and viscoplastic rheology. Compressibility promotes an irregular, more sluggish, flow field and encourages greater interaction between upper and lower boundary instabilities (Steinbach et al., 1989). Viscous heating redistributes buoyancy sources (Jarvis and McKenzie, 1980) and latent heat can reverse phase boundary distortion caused by the advection of ambient temperature (Schubert et al., 1975). These non-Boussinesq effects could play an important role in determining the lateral variations in temperature and phase that may give rise to rapidly varying waveforms. Using constraints from experimental and theoretical mineral physics we predict seismic wavespeed structure from the temperatures and phases determined by the convection calculations. Finally, we compute synthetic seismograms to analyze S, SdS, and ScS arrivals.

Equations
We employ the truncated anelastic liquid approximation (TALA) for infinite-Prandtl-number flow (Jarvis and McKenzie, 1980;Ita and King, 1994) with a divariant phase change using CitcomS (Zhong et al., 2000;. The conservation equations of mass, momentum, and energy (nondimensional) are: where ρ is density, u velocity, P dynamic pressure, τ deviatoric stress tensor, Ra thermal Rayleigh number, α thermal expansivity, T temperature, Rb phase Rayleigh number, Γ phase function, g gravitational accceleration, η viscosity, t time, c heat capacity, κ thermal diffusivity, T S surface temperature, Di = α 0 g 0 R 0 /c 0 dissipation number, and H internal heat production rate. Depth-dependent reference state parameters are denoted by overbars, and the r subscript denotes a unit vector in the radial direction. Phase buoyancy and latent heat are included using a phase function Γ, effective heat capacity c , and effective thermal expansivity α (Richter, 1973;Christensen and Yuen, 1985), defined as: where (p Γ , T Γ ) is the reference pressure and temperature of the phase transition, γ Clapeyron slope, d phase transition width, and p h hydrostatic pressure. The phase function Γ = 0 for Pv and Γ = 1 for pPv.
The thermal Rayleigh number, Ra is defined by dimensional surface values: Our Ra is therefore about an order-of-magnitude larger than a definition with mantle thickness (Table 1).

Rheology
We use a temperature-dependent viscosity, governed by an Arrhenius law, and with plastic yielding to generate a mobile upper surface and strong slabs (Moresi and Solomatov, 1998;Tackley, 2000): The reference viscosity η 0 is defined at the CMB (T = 1). We define the strength envelope of the lithosphere using Byerlee's Law (Byerlee, 1978): where r is the radius, σ yield stress, a cohesion, b brittle stress gradient, and σ y ductile yield stress. The effective viscosity η ef f is computed by: where˙ is the second invariant of the strain rate tensor. The cohesion is set to a small value a = 1 × 10 4 , and remains constant to limit the number of degrees of freedom. We define a yield stress envelope characterized by brittle failure from the surface to 160 km depth and ductile failure for greater depths (b = 40σ y ). Preliminary models reveal b = 1.6×10 9 produces plate-like behavior with strong downwellings.

Model setup
We solve the conservation equations in a 2-D cylindrical section (2 radians) with free-slip and isothermal upper and lower surfaces and zero heat flux sidewalls. Internal heating is neglected (H = 0). Radial mesh refinement provides the highest resolution of 14 km within the thermal boundary layers and the lowest resolution of 40 km in the mid-mantle. Lateral resolution is 0.45 degrees (50 and 28 km at the top and bottom, respectively).
We construct a depth-dependent reference state for density using a polynomial fit to PREM in the lower mantle (Solheim and Peltier, 1990;Dziewonski and Anderson, 1981). The pressure-dependence of the thermal expansion coefficient is derived from experimental data for MgSiO 3 (Mosenfelder et al., 2009) and varies by approximately a factor of 4 across the mantle. Gravity and heat capacity are approximately constant throughout the mantle and therefore we adopt uniform profiles. Thermal conductivity is likely dependent on pressure and temperature, although there are significant uncertainties, particularly in multi-phase and chemically complex systems. We choose to neglect this complexity by assigning the same thermal conductivity to Pv and pPv independent of pressure.
The initial condition is derived from a precalculation (extended-Boussinesq) with similar parameters (Table 1) while excluding the phase transition. This reduces the spin-up time for the TALA model. We integrate the model for several 10,000s of time steps (dimensionally several Gyr) such that the initial condition no longer affects our results.

Post-perovskite phase transition
We set the ambient phase transition absolute temperature (T Γ + T S = 2600 K) to the horizontally averaged temperature 300 km above the CMB in the precalculation. The free parameters are the ambient phase transition pressure p Γ and the Clapeyron slope γ, which we vary to explore a range of dynamic scenarios and consequences for seismic wave propagation. The Clapeyron slope is varied between 7.5 MPa K −1 and 13.3 MPa K −1 , which is consistent with geophysical constraints (Hernlund and Labrosse, 2007) and mineral physics theory and experiment (see review by Shim, 2008). Sidorin et al. (1999b,a) favor a Clapeyron slope of 6 MPa K −1 , although larger values fit the data nearly as well.
The interaction of a simple conductive geotherm with the phase transition can produce a double crossing of the phase boundary (Hernlund et al., 2005). Ambient mantle enters the pPv stability field in D and then transforms back to Pv just above the CMB. In this scenerio a double crossing occurs when the phase transition temperature at the CMB is less than the CMB temperature. This is in part because a conductive temperature profile increases monotonically with pressure. However, temperature profiles that sample thermal heterogeneity in D (slabs, plumes) can intersect the pPv phase space multiple times even for phase transition temperatures that are greater than the CMB temperature. The CMB is an isothermal boundary so we can always compute the stable phase at the base of the mantle for each model a priori (Table 2). If pPv is stable at the CMB, at least one phase transition must occur independent of temperature. If Pv is stable, a double crossing is facilitated for mantle regions in the pPv stability field.

Seismic waveform modeling
We express the relative perturbation to a parameter X using the notation δX = ∂lnX ≡ ∆X/X where ∆X is the change in X. The perturbation to the P-wave velocity, δvp, and S-wave velocity, δvs, can be expressed as : where K is the adiabatic bulk modulus, G shear modulus, ρ density, and R 1 = G/K is equal to 0.45, which is similar to PREM at 2700 km depth. δρ is computed directly from the geodynamic calculations (Appendix A). The perturbations to the elastic moduli are decomposed into phase and thermal parts, assuming they are linearly independent: where dT is a temperature perturbation relative to a reference adiabat. We derive the reference adiabat semi-empirically by considering the horizontally averaged temperature field for the models with a phase transition (Fig. 1). Since a phase transition with a positive Clapeyron slope enhances convection we exclude the model without a phase transition because the adiabatic temperature gradient is less steep and the thermal boundary layers thicker, although the differences are not significant. We construct a reference adiabat with a foot temperature of 0.47 (dimensionally 1732 K) using the dissipation number and thermal expansion profile defined for the geodynamic calculations. These temperatures are contained within the range of time-average geotherms from the models (grey region in Fig. 1) between 500 and 2250 km depth. Temperature perturbations are then determined by removing the reference adiabat from the computed temperatures. We use ∂lnK/∂T = −0.1217 and ∂lnG/∂T = −0.319 derived from a theoretical calculation for MgSiO 3 (Oganov et al., 2001). However, we apply these values to both the Pv and pPv phase. Note that these derivatives are non-dimensional because T is a non-dimensional temperature.
Across the Pv-pPv transition, we prescribe fractional changes in both the shear wave and the compressional wave and compute the corresponding derivatives (Appendix A): where δvs Γ and δvp Γ are fractional perturbations to the S-and P-wave velocity, respectively, due to the pPv phase transition. Beneath the Cocos Plate the average S-wave increase is approximately 2% (δvs Γ = 0.02) and the P-wave variation is a few fractions of a percent (Hutko et al., 2008). For simplicity we assign δvp Γ = 0 but note that short period P-wave reflections are observed in other regions such as beneath northern Siberia (e.g., Weber, 1993). The 1-D seismic model IASP91 (Kennett and Engdahl, 1991) converts velocity perturbations to absolute values. We generate synthetic seismograms using the WKM technique to resolve the SdS triplication originating from high-velocity regions of D (Ni et al., 2000). In this initial study we trace rays through seismic models computed from 1-D temperature and phase profiles extracted from the 2-D convection calculations. We select a Gaussian source time function with 4 s duration and bandpass the synthetic data between 4 and 100 s.

Overview of D slab dynamics
We describe typical slab dynamics using model S2 (Fig. 2) as all models exhibit similar behavior. The pPv phase boundary is 300 km above the CMB and the Clapeyron slope is 7.6 MPa K −1 . Initially, the upper boundary layer thickens and the CMB region warms ( Fig. 2A). The pPv region has constant thickness except where a few thin stationary plumes emanate from the lower thermal boundary layer and the phase boundary is perturbed to higher pressure. A slab forms from the upper thermal boundary layer when the yield stress is exceeded; sometimes the downwelling is generated by several instabilities that unite to form a dominant downwelling. As the slab descends it warms through adiabatic and viscous heating in addition to diffusion.
Pre-existing plumes at the CMB are displaced laterally by the downwelling, such as the plume at the tip of the slab (Fig. 2B). The cold slab geotherm (Fig. 2G, black solid line) intersects the phase boundary (grey dashed line) at lower pressure than ambient mantle, which warps the phase boundary upward producing a thick pPv layer. Latent heat release acts to restore the phase boundary to its unperturbed position and generate a body force that opposes the motion of the slab. However, the amplitude of the thermal anomaly and phase boundary deflection from advected temperature generate the dominant buoyancy forces that are largely unmodified by latent heating.
The slab smothers the CMB, trapping a portion of hot mantle that begins to vigorously convect. During this period the temperature-dependent rheology causes the slab to deform more easily as it warms. The trapped mantle accumulates buoyancy, forming a plume that pushes up the overlying slab (Fig. 2C). As the hot upwelling rises from the CMB it enters the Pv stability field which creates strong phase heterogeneity and provides additional positive buoyancy because it is surrounded by pPv. However, latent heat absorption cools the plume and retards the deflection of the phase boundary, although these effects appear negligible. The upwelling erupts through the slab and suppresses the height of the phase boundary (Fig. 2D). The upwelling loses its thermal buoyancy through adiabatic cooling and thus its temperature anomaly (relative to the reference adiabat) decreases rapidly as it rises.
Secondary plumes form from patches of thickened boundary layer on either side of the original upwelling (Fig. 2D). These also punch through the slab and entrain some of the remaining cooler material to create patches of deformed slab separated by plume conduits (Fig. 2E, but more clearly visible in Fig. 2O). A second downwelling at the rightmost edge of the domain displaces these patches and thermal contrasts create topography on top of the pPv layer. As heat diffuses into the slab the thermal anomaly is eradicated and the pPv layer returns to uniform thickness punctuated by a few plumes as in Fig. 2A. This process is aperiodic.

Overview of seismic wavespeed anomalies
We compute P-and S-wave velocity anomalies using the geodynamic data (Fig. 2) assuming that the P-wave is unperturbed by the phase boundary (δvp Γ =0%) and the S-wave experiences a step-wise increase in velocity (δvs Γ =2%) (Fig. 3). P and S phases with similar paths in the lower mantle will therefore generate different waveform complexity. Furthermore, seismic profiles depend on the evolutionary stage of the slab and its interaction with the thermal boundary layer and pPv phase transition (Fig. 3).
At 0 Myr, the seismic profiles are largely unperturbed except for the pPv phase boundary (Fig. 3P) and reduced basal velocity caused by the lower thermal boundary layer. The slab introduces thermal heterogeneity that produces a strong seismic gradient across its top surface (Fig. 3Q). This thickens the pPv layer and generates a much steeper velocity gradient for the S-wave than just the pPv transition alone (Fig. 3Q). In this region the wavespeed reduction above the CMB is also larger than it would be without an overlying slab.
The growth of a plume beneath the slab produces four distinctive seismic velocity gradients (Figs. 3H,R). Across the top of the slab the gradient is positive and steep, but wavespeeds rapidly decline in the plume due to temperature and phase (S-wave only). Wavespeeds are consistently low in the plume with a gradient comparable to the 1-D model. The thermal boundary layer at the base of the mantle generates a second negative velocity gradient, although the actual decrease in velocity is significantly less than the transition from slab to plume.
When the plume punches through the slab the rising hot material suppresses wavespeeds and generates seismic profiles similar to those at 0 Myr (compare Figs. 3I,S with Figs. 3F,P). The plume head is surrounded by cooler remnant slab which introduces a negative wavespeed gradient at 2100 km depth. Secondary plumes develop in the lower thermal boundary layer and break apart the remnant slab. Velocity profiles through a developing plume in the pPv stability field (Figs. 3J,T, black solid line) exhibit a low negative gradient at the CMB. This contrasts to a Pv plume otherwise within a pPv layer which produces a positive velocity gradient except for the very base of the mantle (compare Figs. 3J,T with Figs. 3H,R, black solid line). The neighboring blocks of cool deformed slab produce a sharp phase boundary with relatively high wavespeeds (Figs. 3J,T, black dashed line). The lower thermal boundary is suppressed beneath the slab, which generates a high negative velocity gradient at the base.

Dynamic and seismic variations
Our models explore variations in the phase boundary height and Clapeyron slope. Increasing the pPv transition pressure at fixed temperature is equivalent to reducing the pPv transition temperature at a fixed pressure because the pPv stability field is simply translated in P-T space.
The leftmost columns of Figs. 4 and 5 show models with an increased Clapeyron slope of 10.6 MPa K −1 . Models S3, S4, and S5 have equivalent transition pressures of 300 km, 450 km, and 150 km above the CMB, respectively. We show the S-wave velocity but the perturbations to the P-wave are comparable except for the step-wise increase due to the pPv phase transition. Model S3 (Figs. 4A,E,B,F) has a transition pressure of 300 km above the CMB which is the same as the "typical" model in Section 4.1. However, the Clapeyron slope is larger so the stable phase at the CMB is now Pv rather than pPv and plumes form at the CMB in the Pv stability field (Figs. 4A,E, black solid line). This accentuates the wavespeed contrast with the overlying cooler pPv by producing a lower velocity basal layer (compare Figs. 3O,T with Figs. 5A,E, black solid line).
Figs. 4A,E reveal lower boundary layer instabilities at various stages of their life cycle. This includes thickening of the boundary layer, plume detachment from the CMB, and eruption through the pPv layer with a trailing conduit. Plume heads become Pv within the pPv layer. Seismic waves may therefore sample all of this complexity within small ranges of azimuth or epicentral distance. At later time, several plumes punch through the slab and generate separate blocks of pPv that are surrounded at the base and edges by hot Pv (Figs. 4B,F).
Model S4 has a lower transition pressure than S3 (450 km versus 300 km) which is equivalent to translating the phase boundary to higher temperature (Figs. 4C,G) such that the stable phase at the CMB is pPv. Plumes form in the pPv stability field and intersect the Pv-pPv phase boundary as they rise through D . This forms isolated pockets of hot Pv within the pPv-dominated layer (Fig. 4G, black dashed line) which causes multiple perturbations to the S-wave velocity (Fig. 5G, black dashed line). The thermal profile crosses the phase boundary multiple times and is an example of a "triple crossing" (Fig. 4G, black dashed line). Escaping plume heads also create concave-up notches in the phase boundary at the top of D (Figs. 4C,G, black solid line). This is also observed in S2 (Figs. 2O,T) but is more noticeable here because the Clapeyron slope is larger. The pPv layer has a smooth upper surface at longer wavelength, which contrasts with models that produce distinct pPv blocks. However, the layer is permeated by plumes that produce low-velocity regions extending radially from the CMB.
Model S5 has a higher transition pressure than S3 (150 km versus 300 km above the CMB) but behavior is comparable because both models facilitate a double crossing for the average 1-D geotherm. However, perturbations to the phase boundary height are more comparable to the ambient thickness of the pPv layer (Fig. 4D). Blocks of pPv are more separate and there is relatively larger topography on their tops.
A larger Clapeyron slope (13.3 MPa K −1 ) enables greater lateral phase heterogeneity (Figs. 4I,J,K). In model S8, pPv only forms within cold slabs in D because the average 1-D geotherm is too warm to sustain a permanent pPv layer (Fig. 4I,M, black solid line). For example, pPv occurs around the tip of a slab as it penetrates D (Fig. 4I) causing substantial variations in phase boundary height (Figs. 4M,N). Plumes can penetrate through older slabs whilst downwelling material continues to accumulate nearby, which forms both regions of pPv blocks and a continuous layer (Fig. 4K).
Finally, a folded slab at the CMB has considerable complexity. Two negative temperature perturbations are generated by the uppermost and lowermost parts of the slab fold with ambient material sandwiched between (Fig. 4L). The top of the lowermost thermal anomaly accompanies the transition to pPv. These features are clearly evident in the seismic profile (Fig. 5P, black dashed line). Seismic velocity gradients are steeper in the cold regions, with the second perturbation accentuated by the phase transition. The negative velocity gradient at the base of the mantle is large because of the overlying slab and Pv is the stable phase at the CMB.

Waveform predictions
Thermal and phase heterogeneity produce wavespeed perturbations that can be explored through waveform modeling. We produce 1-D synthetic seismograms for the S-wave velocity profiles in Figs. 3 and 5. First we overview basic features that are necessary to interpret the synthetic data.
The seismograms are aligned on the IASP91-predicted S arrival and the triplicated SdS phase arrives before ScS. IASP91 (grey seismograms) does not produce a SdS arrival as it lacks a velocity discontinuity just above the CMB. The height of the interface can be determined by SdS-S and ScS-SdS differential travel times. For a given epicentral distance, SdS-S increases and ScS-SdS decreases as the phase boundary is perturbed to high pressure. The relative amplitude response of SdS depends on the magnitude of the velocity increase across the discontinuity and ScS-S constrains the volumetric velocity anomaly within D (below the interface). SdS overtakes S as the first arrival at post-crossover distances, which is ∼ 85 degrees for the 1-D seismic model SLHO for a focal depth of 600 km (Lay and Helmberger, 1983). A negative seismic velocity gradient above the CMB slows SdS which pushes the postcrossover SdS and S arrivals together. Furthermore, more energy turns below the discontinuity, which weakens the post-crossover S phase (Young and Lay, 1987).
A recently subducted slab produces a large velocity increase at the phase boundary because of the coincident strong thermal gradient and phase jump (Fig. 3L) which generates large SdS amplitude (Fig. 6B). Conversely, for a slab that has warmed at the CMB the thermally induced velocity jump at the D discontinuity is less and the SdS amplitude is reduced (Fig. 6E). Relative to a thinner D region, a slab with a broad thermal anomaly elevates the phase boundary, which causes SdS to appear at shorter distances, decreases SdS-S, and reduces the SdS-S crossover distance (compare Fig. 6B with Fig. 6E). Thick D regions with cold slabs have a large positive volumetric wavespeed anomaly which causes ScS to arrive early (ScS-S decreases) and therefore reduces the ScS-S crossover distance (Fig. 6B). As D becomes increasingly thin, ScS arrives closer to the IASP91 prediction, SdS-S increases, and ScS-SdS decreases until SdS arrives as a precusory shoulder to ScS (Figs. 8D,E). Post-crossover SdS amplitude is large for fast regions that extend from the phase boundary to the CMB (e.g., Fig. 6B).
A thickening boundary layer at the base of a slab has a negative wavespeed anomaly which increases the arrival time of ScS (ScS-S increases) by reducing the volumetric wavespeed anomaly (Fig. 6F). As a plume develops, the wavespeeds beneath the D discontinuity continue to decrease and ScS arrives with the IASP91 prediction (Fig. 7G) or is delayed ( Fig. 6C and Fig. 7B). These results are consistent with the volumetric velocity anomaly controlling ScS-S travel times (e.g., Sidorin et al., 1999a). The temperature anomaly at the phase interface decreases as the plume heats the overlying slab, which reduces SdS amplitude by decreasing the velocity jump (Fig. 7C). However, SdS-S does not vary because this is sensitive to the height of the discontinuity above the CMB, which is unaffected by plume growth within D until the plume erupts through the phase boundary. A warm plume bends rays away from the normal which reduces S and SdS amplitude beyond the SdS-S crossover relative to a horizontal slab (compare Fig. 7G and Fig. 7H beyond ∼ 80 degrees). Fig. 7F shows the seismogram for a Pv plume head within a pPv layer (Fig. 4C, black dashed line). SdS arrives from the top of the pPv layer but a second triplicated arrival from the warm plume to pPv transition at higher pressure is absent. The slow velocity within the plume presumably creates a shadow zone that suppresses this second arrival. The waveforms for this profile are similar to those produced by a simple phase transition overlying a thermal boundary layer (Fig. 7A). Although the velocity structure differs from IASP91, the arrival time of S and ScS is similar.
The seismogram for a plume erupting through a pPv layer has an additional (weak) triplication phase that arrives after S for distances 60-68 degrees (Fig. 6D). This phase originates from cooler slab material carried above the plume head. Beyond 75 degrees the direct S arrival enters a shadow zone due to ray bending from sampling the slower velocity plume. Since the wavespeed anomalies are generally negative throughout D , the SdS-S and ScS-S crossover distances are > 84 degrees. Two triplicated arrivals are also produced by a folded slab interacting with the phase transition (Figs. 4L,5L). The first arrives between 62-70 degrees and originates from the top of the uppermost slab fold (Fig. 8H). The second arrives beyond 65 degrees and is SdS from the velocity increase caused by the slab interacting with the phase boundary. Combined analysis of multiple triplicated arrivals, possibly with other phases, may be able to determine the characteristic lengthscale of slab folding above the CMB.
A simple thermal boundary layer at the CMB delays ScS, and marginally delays the arrival time for S at large distances (> 80 degrees) (Figs. 8A,F). A warm conduit emanating from the CMB delays S beginning at shorter distances (> 76 degrees) and further suppresses the amplitude of the arrival (Fig. 6A). Additionally, a weak SdS phase is faintly visible beyond about 72 degrees and is generated by the small velocity increase at the phase transition (Fig. 6A). ScS-SdS is small because the discontinuity is close to the CMB.

Discussion
Slabs rapidly advect and blanket the CMB from episodic flushing of the upper thermal boundary layer. The maximum thermal anomaly of slabs at the CMB is 0.33 (1200 K), which is at the upper bound of estimates (Steinbach and Yuen, 1994). Additionally, the simple two-phase model (Pv and pPv) and uncertainties in material properties preclude an accurate mapping of temperature and phase to seismic wavespeed. However, we focus on changes in seismic velocity gradient because the absolute magnitude of velocity perturbation is not critical for analyzing trends in seismic triplication data. Similarly, choosing an alternative reference geotherm would not significantly affect our results. The largest discrepancy between the horizontally averaged temperature and reference adiabat is in the region of subadiabatic temperature gradients above the CMB. If we choose a horizontally averaged temperature as reference (extrapolating across the boundary layers) the wavespeed anomaly for slabs (plumes) will slightly decrease (increase) in amplitude. Finally, diffuse thermal anomalies at the CMB would not produce sharp seismic gradients and lateral variation in the phase boundary necessary to investigate high-frequency waveform complexity.
Downwellings displace pre-existing plumes at the CMB because low-viscosity conduits are established between subduction cycles. Therefore plumes infrequently develop at slab tips in comparison to models with an unperturbed lower thermal boundary layer Tan et al., 2002). Episodic subduction promotes plume development beneath slabs because existing basal instabilities are swept aside. Slabs insulate the CMB allowing the boundary layer to thicken and form upwellings. Tan et al. (2002) investigate these dynamics and propose two types of outcomes: "normal plumes" and "mega-plumes". Qualitatively the plumes that we observe are more akin to "normal plumes" because they are not voluminous eruptions and dissipate quickly above the CMB. We anticipate different behavior because our models are compressible with a phase transition, whereas the fine-scale models of Tan et al. (2002) are Boussinesq. Convective vigor and adiabatic cooling control the temperature anomaly of plumes in compressible flow models (Zhao and Yuen, 1987). After plumes dissipate the thermal anomaly of the slab they combine into a few single upwellings (Vincent and Yuen, 1988).
Nevertheless, the seismic profile for a mega-plume beneath an ancient slab (Fig. 11, right panel in Tan et al., 2002) is comparable to a plume within a pPv layer (e.g., Fig. 3M,R). The later includes an S-wave velocity increase across the Pv-pPv phase boundary. Uncertainties when mapping temperature and phase to seismic wavespeed prevent 1-D seismic analyses from identifying a warm plume (seismically slow) versus a hot mega-plume (slower). Furthermore, we extract 1-D temperature and phase profiles from 2-D convection calculations and are therefore insensitive to the total volume of the plume. The profiles above the CMB are different because Tan et al. (2002) compute δvs using the horizontally averaged temperature field as a reference geotherm, whereas we use an adiabat. Our models therefore include a negative seismic velocity gradient just above the CMB in accordance with seismic data (Cleary, 1974).
The models suggest several mechanisms that may explain the origin of some seismic scatterers. First are small-scale plumes at different stages of development forming in a pPv layer (Fig. 4A). The phase difference between a Pv plume and the surrounding pPv layer will further increase wavespeed variations in comparison to just the temperature effect. Second are steepsided isolated blocks of pPv (Fig. 4B). Third are concave-up notches in the phase boundary caused by escaping plume heads (Fig. 4C). Furthermore, these processes often operate in combination.
The synthetic seismograms reveal triplicated arrivals originating from the D discontinuity and large thermal gradients, e.g., a folded slab (Fig. 8H). The timing and amplitude response of the arrivals provide insight into the velocity increase, height of the phase transition above the CMB, and approximate thickness of the fast region. However, S, SdS and ScS in the distance range 60-83 degrees are largely insensitive to small-scale features beneath the phase boundary. Negative seismic gradients are inherently difficult to resolve because energy is not refracted to the surface, which means that low-amplitude reflected arrivals, often below the noise level, have to be analyzed (e.g., Flores and Lay, 2005). However, stacking techniques can be used to enhance low-amplitude arrivals (e.g., Kito and Krüger, 2001;Kito et al., 2007). The ScS-S differential travel time constrains the volumetric wavespeed anomaly of D but is not sensitive to the distribution of heterogeneity.
Seismic anisotropy is expected to modify our analysis. Notably, the magnitude and sign of the seismic velocity perturbation depends on the wave propagation direction for aligned pPv (e.g., Ammann et al., 2010;Thomas et al., 2011). This can modulate the amplitude of seismic phases. However, since our geodynamic models do not track strain accumulation and the development of fabric (e.g., McNamara et al., 2002) we do not consider anisotropic effects in this study. Furthermore, future models can explore lateral variation in wavespeed and the focusing and defocusing of seismic rays due to topography on the phase boundary.

Conclusions
The models reveal complex interaction of slabs, plumes, and the Pv-pPv transition which produces significant thermal and phase heterogeneity in D over small distances. Slabs deflect the phase boundary and disturb the lower thermal boundary layer, pushing aside pre-existing upwellings. Plumes regularly develop beneath slabs and distort the phase boundary and slab morphology as they erupt from the CMB. Plume formation occurs less frequently at slab edges because a few stationary upwellings drain the boundary layer between cycles of subduction. Advected isotherms produce more complicated thermal structure than conductive profiles. This enables multiple crossings of the phase boundary, e.g., a Pv plume head within a pPv layer can produce a "triple crossing".
The topology of the pPv layer depends on the transition pressure (or equivalently transition temperature since both translate the phase boundary in P-T space) and the Clapeyron slope. The models reveal a variety of behavior including: (1) distinct pPv blocks separated by upwellings, (2) notches at the top of a pPv layer caused by plume heads, (3) regions of Pv embedded within a pPv layer due to upwellings. Furthermore, these scenerios have interesting consequences for seismic scattering phenomena which will be explored in a future study.
In synthetic seismograms, we cannot resolve small-scale complexity beneath the D discontinuity using solely S, SdS, and ScS triplication data to 84 degrees. This is because of inherent limitations in resolving negative seismic velocity gradients. In future work we will model the seismic wavefield using a finite difference approach to explore SH-, SV-, and P-waves for many D arrivals (e.g., S, Scd, ScS, SKS, SKKS) and potentially reverberations in the pPv layer. Extending the distance range to 110 degrees will also allow us to analyze diffracted phases (e.g., Sd) and arrivals beyond the SdS-S crossover distance.
There are possibilities for developing more complex dynamic models. Compositionally stratified slabs, possibly interacting with a chemical layer at the CMB, will also generate waveform complexity (Tackley, 2011). Furthermore, chemical heterogeneity in D influences seismic velocities and the Pv-pPv transition pressure and depth of the mixed phase region. The effects of anisotropy need evaluating. Thermal conductivity strongly controls heat transport at the CMB and therefore the thermal evolution of slabs.
M. Jackson for comments. We also thank Allen K. McNamara and an anonymous reviewer for constructive comments. This work was supported by NSF grant EAR-0855815, and DS was supported by a CIW/DTM postdoctoral fellowship. Figures were produced using GMT (Wessel and Smith, 1998).   Figure  S1 No phase transition S2 0.403 (300) 0.114 (7.6) pPv 2, 3, 6 S3 0.403 (300)   The pPv layer at the CMB is contoured in black and the phase transition parameters (height above the CMB and Clapeyron slope) are given below each panel. Black solid and dashed lines represent profile locations. (E-H) and (M-P) correspond to differential temperature profiles at the profile locations for A-D and I-L, respectively. The grey solid line denotes the reference adiabat (defined zero) and the grey dashed line shows the Pv-pPv phase boundary. Given times are relative to 0 Myr for that particular model.   where δvs Γ and δvp Γ are the fractional perturbations to the S-and P-wave velocity, respectively, due to the Pv-pPv phase transition. By substitution: (A.14)