A geodynamic and mineral physics model of a solid-state ultralow-velocity zone

Recent results (Wicks et al., 2010) suggest that a mixture of iron-enriched (Mg,Fe)O and ambient mantle is consistent with wavespeed reductions and density increases inferred for ultralow-velocity zones (ULVZs). We explore this hypothesis by simulating convection to deduce the stability and morphology of such chemically-distinct structures. The buoyancy number, or chemical density anomaly, largely dictates ULVZ shape, and the prescribed initial thickness (proxy for volume) of the chemically-distinct layer controls its size. We synthesize our dynamic results with a Voigt-Reuss-Hill mixing model to provide insight into the inherent seismic tradeoff between ULVZ thickness and wavespeed reduction. Seismic data are compatible with a solidstate origin for ULVZs, and a suite of these structures may scatter seismic energy to produce broadband PKP precursors.


Introduction
The large chemical, density, and dynamical contrasts associated with the juxtaposition of liquid iron-dominant alloy and solid silicates at the coremantle boundary (CMB) are associated with a rich range of complex seismological features. Seismic heterogeneity at this boundary includes small patches of anomalously low sound velocities, called ultralow-velocity zones (ULVZs). Their small size (5 to 40 km thick) (e.g., Garnero and Helmberger , 1996) and depth (>2800 km) present unique challenges for seismic characterization.
ULVZs were first noted with teleseismic SPdKS phase diffracting along the CMB to determine the P-wave velocity at the base of the mantle beneath the central Pacific (e.g., Garnero et al., 1993;Garnero andHelmberger , 1995, 1996;Helmberger et al., 1996). The compressional wavespeed decreases by 5 to 10% over a depth of 5 to 40 km above the CMB, consistent with PcP precursors (e.g., Mori and Helmberger , 1995;Revenaugh and Meyer , 1997;Hutko et al., 2009). In other locations, ULVZs are seismically absent or below the detection threshold, such as the North Pacific (Rost et al., 2010b), which may imply they are not globally ubiquitous (see Thorne and Garnero, 2004). ULVZs have been correlated with the location of hotspots (Williams et al., 1998) and tentatively to the edges of large low shear velocity provinces (e.g., Lay et al., 2006). Recent thermochemical convection calculations lend support to these spatial correlations (McNamara et al., 2010).
Precursors and postcursors to the converted core-reflected phase ScP can potentially constrain the P and S wavespeed, thickness, and density of ULVZs (e.g., Garnero and Vidale, 1999;Reasoner and Revenaugh, 2000). Based on this approach, analyses using small-aperature and short-period arrays have elucidated the structure between Tonga-Fiji and Australia at high resolution Idehara et al., 2007;. These studies report P and S wavespeed reductions of ≈ 8% and ≈ 24% respectively, a thickness of about 10 km, and a density increase of around 10%.
All of these studies characterize structure by a 1-D model. Although quasi-1D models utilize different structures for the source and receiver paths (e.g., Garnero and Helmberger , 1996;Helmberger et al., 1996), only a few 2-D models have been reported. Probing the CMB beneath the southwest Pacific, Wen and Helmberger (1998b) model SKS-SPdKS observations with Gaussian-shaped ULVZs of approximately 40 km in height, 250 to 400 km lateral extent, and a P-wavespeed drop of about 10%. PKP precursors suggest that these larger structures are composed of smaller undulations (Wen and Helmberger , 1998a). A concave-down upper interface is necessary to model structures beneath Africa and the eastern Atlantic . In both 1-D and 2-D models there are tradeoffs between the height and velocity decrease of ULVZs (e.g., Garnero and Helmberger , 1998;Wen and Helmberger , 1998b).
A partial melt origin for the ULVZs predicts a P to S wavespeed reduction of about 1:3 (Williams and Garnero, 1996;Berryman, 2000;Hier-Majumder , 2008). In this model, the velocities of the average assemblage are decreased by melt formed either by fluid reaction products of Fe liquid with silicate mantle or by partial melting of Fe-rich mantle. The hypothesis is consistent with the correlation between ULVZs and hot spots (Williams et al., 1998) and the broad agreement with P to S velocity reductions determined from core-reflected phases. Early dynamical calculations question the ability to produce a dense and non-percolating melt phase in the deep Earth (Hernlund and Tackley, 2007). However, the stirring of ULVZs by the larger-scale convective motions of the mantle can potentially maintain a partially molten region (Hernlund and Jellinek , 2010). Partial melt may also be trapped within ULVZs because of textural changes (e.g., , but theoretical and experimental justifications are incomplete, particularly at the temperatures and pressures of the CMB region. Iron enrichment of solid phases, specifically the increase in Fe/(Fe+Mg) ratio, can simultaneously increase density and reduce compressional and shear velocity (e.g., Karato and Karki, 2001). This partly inspired the notion of solid, iron-rich ULVZs, such as a metal-bearing layer (Knittle and Jeanloz , 1991;Manga and Jeanloz , 1996), subducted banded iron formations (Dobson and Brodholt, 2005), or iron-enriched post-perovskite (Mao et al., 2006;Stackhouse and Brodholt, 2008). Iron-rich systems are typically denser than the surrounding mantle, which is required to explain the locations of ULVZs at the base of the mantle.
Recent results show that the sound velocities of a solid iron-enriched (Mg,Fe)O at CMB pressures are low enough, such that only a volumetrically small amount (∼10-20%) mixed with coexisting silicates are needed to explain ULVZs (Wicks et al., 2010). We are motivated to generate a numerical convection model of a thin, iron-enriched (Mg,Fe)O-containing layer interacting with the lowermost thermal boundary layer. Mantle convection thickens (thins) a dense layer beneath upwellings (downwellings) and controls its spatial distribution (Davies and Gurnis, 1986). A persistent stable layer has a density contrast of a few percent (e.g., Garnero and McNamara, 2008) and may convect internally (Hansen and Yuen, 1988). We explore whether these characteristics are manifested in this iron-enriched (Mg,Fe)Ocontaining layer. We determine the steady-state morphology of such a layer and develop an integrated and self-consistent ULVZ model which uses constraints from geodynamics and mineral physics and is consistent with seismic data.

Equations and solution methods
We apply the Boussinesq approximation to model thermochemical convection using CitcomS (Zhong et al., 2000;Tan et al., 2007) because the pressure range of our domain is small and energy dissipation is negligible with a greatly reduced thermal expansion coefficient at high pressure. The equation for the conservation of mass for an incompressible fluid is: where u is velocity. The non-dimensional momentum equation is: where P is the dynamic pressure, η is the viscosity,˙ is the deviatoric strain rate, Ra is the thermal Rayleigh number, T is temperature (nondimensionally ranging from 0 to 1), Rb is the chemical Rayleigh number, and r is the radial unit vector. Composition, C, ranging from 0 to 1, encompasses two chemical components. The first component (C = 0) is ambient mantle, the second (C = 1) is a two-phase mix containing iron-enriched (Mg,Fe)O that constitutes the distinct chemistry of the ULVZ. We equivalently refer to the latter component as the chemical, or (Mg,Fe)O-containing, component. The Rayleigh numbers are defined as: where the parameters and their values are provided in Table 1. A useful non-dimensional parameter is the buoyancy number, B, which describes the chemical density anomaly normalized by the maximum thermal density anomaly, or equivalently the ratio of chemical to thermal buoyancy: The quantity of heat-producing elements in the lower mantle and the dynamic implications for small-scale structures such as ULVZs remain uncertain. For simplicity, we solve the non-dimensional energy equation for temperature without an internal heat source: where t is time.
The equation for chemical advection is: We represent the ULVZ component by a set of tracer particles (initially about 100/cell) that are advected using a predictor-corrector scheme (McNamara and Zhong, 2004). The volume fraction of ULVZ material is proportional to the absolute local concentration of tracers with a truncation applied to prevent unphysically large values (see Tackley and King, 2003, the truncated absolute method). In comparison to the ratio method (Tackley and King, 2003), this technique provides greater computational speed because fewer total tracers are required. Although the number of tracers per cell is greater for the absolute method to achieve the same resolution, a large volume fraction of ambient material in our models is devoid of tracers. By contrast, the ratio method requires tracers to fill the space of the entire domain.

Domain and rheology
We solve for the flow within a small region near the CMB in a cylindrical geometry. Our domain (513×145 nodes) spans approximately 23 • with a height of 488 kilometers above the CMB. The radial resolution is 2 km in the lowermost 256 km of the domain. Above 256 km the grid spacing increases incrementally by 2 km to a maximum of 20 km.
Non-dimensional viscosity (normalized by η 0 , see Table 1) follows a linearized Arrhenius law; recalling that T is non-dimensional temperature: where q specifies the order of magnitude of viscosity variation within the domain (see Table 1) and is largely restricted by numerical considerations. The stability of the lower thermal boundary layer is strongly dependent on q (e.g., Yuen and Peltier , 1980;Loper and Eltayeb, 1986). The background viscosity modulates the critical size that an instability must attain before a diapir can detach from the thermal boundary layer (e.g., Olson et al., 1987). Therefore, the reference viscosity η 0 (see Table 1) is representative of the unperturbed ambient mantle beyond the influence of the thermal boundary layer.

Boundary and initial conditions
The comparatively large lateral extent of the domain limits the sensitivity of the solution to the insulating sidewalls. The bottom boundary (CMB) is isothermal and free slip. We apply two different top boundary conditions. An isothermal and free slip condition (impermeable) excessively constrains flow and retards upwellings. In contrast, a boundary prescribed by zero normal stress, zero horizontal velocities, and zero heat flux (permeable) (e.g., Tan et al., 2002) lessens resistance so that instability growth accelerates. The impermeable and permeable conditions effectively bound the characteristic dynamic timescales that we expect for an evolving thin boundary at the base of a thick domain (Olson et al., 1987).
Ambient deep mantle is prescribed at zero (non-dimensional) temperature, and the isothermal lower boundary is set to unity. The temperature scaling is given in Table 1. We apply a half-space cooling model with an age of 25 Myr to generate a lower thermal boundary layer with a thickness of ≈66 km. We place a thin layer of the chemical component at the base within the thermal boundary layer. Both the buoyancy number (B) and initial thickness (d ch ) are free parameters (see Table 1). We explore cases where d ch ranges from 32 km to 2 km, which is between 2 and 33 times smaller than the initial thermal thickness. This is reasonable given the expected thickness (5-40 km) of ULVZs.
We allow the chemical layer to be swept into one coherent structure by applying a half sinusoidal perturbation to the temperature field (initially of non-dimensional magnitude of 0.05) both radially and laterally. This permits the outcome of the different numerical experiments to be compared more readily. The size of our domain is carefully specified to focus on individual ULVZ structures and is not designed (nor suitable) to analyze the spatial distribution of many piles. Physically, our chosen initial perturbation is akin to an upwelling on the CMB caused by the (symmetric) downwelling of ambient material on either side.

Transient period
We refer to the initial stage of the model evolution as the transient period. Although we are primarily concerned with the long-term stability and morphology of ULVZs, during this time we make several observations. We overview the evolution of four extreme cases (all impermeable), given by the For (B = 0.5, d ch = 16 km) ( Fig. 1), the thermal boundary layer develops two dominant instabilities within the chemical layer ( Fig. 1A) that are ultimately overwhelmed by the imposed long-wavelength perturbation (Fig. 1B). During this period, the chemical layer develops significant relief caused by the developing plumes, until the plumes merge and a large volume fraction of the layer is expelled from the CMB region (Fig. 1C). The buoyancy of the plume generates viscous stresses that support the relief of the remaining chemical material anchored at the CMB. Stresses reduce as the plume ascends away from the lower boundary and the convection approaches a steady-state, causing a marginal reduction in relief.
The evolution for the case with the same buoyancy number, but a thinner layer (B = 0.5, d ch = 8 km) is similar, but the dynamic timescale is reduced by several tens of million years. This occurs for all cases with a thinner layer, because less thermal buoyancy (less thickening of the boundary layer) is necessary to overcome the intrinsic density increase of the chemical layer.
For the cases with a large chemical density anomaly, (B = 4, d ch = 16 km) and (B = 4, d ch = 8 km), a single plume forms that marginally deflects the layer upward as it leaves the CMB region. The influence of the buoyant plume is lessened as it approaches the top of the domain, which causes the layer to flatten slightly. The large chemical density anomaly prevents substantial topography from forming on the interface.

Steady-state
After the initial transient period, the low effective Rayleigh number ensures that almost all models reach a steady-state, defined by <1 km fluctuations (standard deviation) in structure height. We can thus determine several time-independent characteristics. We record the layer dimensions for the impermeable cases during a steady-state time window defined after the plume reaches the top of the domain but before the downwellings exert sufficient stress to further deform the chemical layer. The permeable cases are unaffected by return flow because material is allowed to escape from the domain.
Five impermeable cases are time-dependent because of an oscillatory flow driven by a pulsating plume that marginally raises and flattens the chemical structure. We compute the standard deviation of the maximum height (km), h sd = 1.5 km. These models are characterized by a low buoyancy number and low initial layer thickness such that the buoyancy of the chemical layer is insufficient to stabilize the flow. We time-average these cases to deduce the stationary steady-state behavior. Fig. 2 depicts all of our parameter space. The relief is defined as the difference between the maximum and minimum height of the structures above the CMB over the lengthscale of the domain (∼1400 km) at steady-state. This metric becomes less appropriate for cases where the ULVZ spreads over a distance greater than 1400 km, thus creating a ubiquitous layer with a much reduced relief even though a potentially voluminous layer exists. Although this may not occur within a larger model domain, a thick continuous layer does not satisfy the small-scale localized nature of the ULVZs. The two largest initial chemical layer thicknesses, d ch = 32 km and d ch = 24 km, generally produce thick ubiquitous layers, or structures with reliefs that are too large to reconcile with seismic observations (i.e. 100 km). We therefore do not discuss these results further.
A certain volume of residual material is necessary to generate structural reliefs within the seismic estimates of 5-40 km. In these models, this can be achieved by either entraining a large amount of material from a thicker chemical layer e.g. (B = 0.5, d ch = 16 km), or by starting with less material e.g. (B = 0.75, d ch = 4 km).
We plot the steady-state morphologies for the impermeable and permeable cases in Figs. 3, 4, and 5. Although the dynamic timescales are affected by the top boundary condition, the long-term evolution of the structures is comparable. Models with low B exhibit triangular morphologies and these shapes progressively flatten as B is increased, with the upper interface becoming more rounded. The morphologies are relatively insensitive to d ch , although for reduced d ch the structures exhibit slightly flatter tops (e.g., Fig. 3) Entrainment modulates the quantity of material available to form the structures for B 0.75. For d ch ≥ 8 km a large volume fraction of the chemical layer is entrained when B = 0.5, which explains the substantially reduced size of the triangular ULVZs. This effect is accentuated for the impermeable cases because of increased viscous stresses. Entrainment is most effective at early model time because viscosities and density differences are comparable (e.g., Sleep, 1988). For this reason, the thinnest initial layer is less affected because it is deeply embedded in the hottest, least viscous part of the domain. Lin and van Keken (2006) describe the entrainment as endmember regimes, which can conveniently be applied to our models. In their regime I, the plume head is formed by ambient material above the CMB, and no chemically-distinct material is entrained. This occurs for B 0.75 in our models. B 0.75 is characteristic of their regime III, particularly for d ch ≥ 8 km, because the dense component predominantly accumulates at the base of the plume with some entrainment by the plume head. In addition to this complex entrainment behavior, relief varies most greatly around B = 1 because of the comparable thermal and chemical buoyancy (Fig. 3).
The half-widths of the structures vary in proportion to B (Fig. 4). Since d ch limits the total volume of material available to form the ULVZ, the halfwidth also varies in proportion to d ch . For B 0.75 the permeable cases produce structures that are slightly less flattened than the impermeable cases, with a subsequent reduction in half-width. The largest variations for both half-widths and relief occur for the thickest initial layer, d ch = 16 km. In Fig. 5 we show that the aspect ratio (relief/half-width) of the structures is strongly dependent on B and only marginally dependent on d ch .

Discussion
PKP precursors provide evidence for scatterers of varying lengthscales at the CMB (Cleary and Haddon, 1972). Our simulations generate structures of various shapes and sizes (e.g., Fig 3) which are comparable to the small-scale and large-scale features used to model broadband data (Wen and Helmberger , 1998a). This suggests that solid-state ULVZs are strong candidates for the origin of the scatterers, particularly if dynamic processes sweep together several ULVZs of differing lengthscales. Many of our dynamic structures also have concave-down upper surfaces that are akin to the 2-D Gaussian-shapes used to model seismic data (e.g., Wen and Helmberger , 1998a,b). A curved surface generates long-period PKP precursors by wide-angle reflection (Wen and Helmberger , 1998a), reduces the amplitudes (relative to a 1-D model) of SKS-SPdKS phases, and produces strong precursors to ScP and ScS (Wen and Helmberger , 1998b). The sharp tops of ULVZs inferred seismically (e.g.,  motivated our a priori choice to define a pre-existing sharp boundary separating the chemical component from ambient mantle (e.g., Fig 6). If the iron-rich (Mg,Fe)O-containing component is formed through a diffusive process, for example by reactions with the core, the interface is not as sharp (e.g., Kellogg and King, 1993).
Countercirculation is observed in all models, and two cases are shown in Fig. 6. Viscous coupling is energetically favorable because it minimizes viscous shear at the material interface (e.g., Richter and McKenzie, 1981). Similarly, in a series of laboratory tank experiments, Jellinek and Manga (2002) observe a preference for material to rise along the sloping interface between a low viscosity chemical layer and ambient material.
In Fig. 6A the ULVZ has sufficient thickness to greatly influence the thermal structure in the lowermost 90 km, and the two-layer convective regime is clearly revealed by Profile II in Fig. 6C. Isotherms are deflected to lower pressure about the central upwelling, which will determine the stable phase assemblages and melt generation zones. Our solid-state model does not require partial melt to account for the dynamic or wavespeed character of ULVZs, but neither is it explicitly excluded. In this regard, partial melt may explain local variations in ULVZ structure, for example beneath the Philippine-Kalimantan region (Idehara et al., 2007).
For B 2, the amplitude of the steady-state perturbation of the chemical layer is much reduced, approximately 35 km in Fig. 6B. The velocity arrows indicate that this flow is generally sluggish with a recirculation of material accommodated by a low viscosity channel at the very base of the domain. Conduction is clearly the dominant heat transfer mechanism (Fig. 6D) and the extreme outcome for this regime is a ULVZ that is subcritical to convection. In this case the ULVZ would act like a 'hot plate' at the base of the mantle and simply offset the thermal boundary layer to marginally lower pressure.
Youngs and Houseman (2009) determine that half-width ∝ B 0.38 and height ∝ B −0.39 based on calculations in a 3-D Cartesian geometry using an isoviscous fluid. In our models, however, we have a temperature-dependent viscosity that introduces a non-linearity which increases the complexity of the convection and questions the applicability of simple scaling relations that involve B. Given this caveat, we find that the exponent of B for half-width ranges between 0.76 and 0.37, and for the height ranges between −0.72 and −0.34.
We now consider the uncertainties in our geodynamic modeling and their implications. The Rayleigh number affects the morphology of a chemical layer and its topography scales approximately in proportion to Ra −1/3 (Tackley, 1998). We assume a thermal conductivity of around 6.6 W/m-K (e.g., de Koker , 2010) but larger values (around 28 W/m-K) have been reported (Hofmeister , 2008). The scaling relation thus predicts that the topography of the chemical layer could be 1.6 times larger than our model predicts. Similarly, an upper bound for lower mantle viscosity is two orders of magnitude larger than our reference value (e.g., Ammann et al., 2010), which by itself could increase relief by a factor of 4.6. Both high thermal conductivity and lower mantle viscosity could collectively elevate the topography by 7.5 times, but mid-range estimates predict a more modest increase of a factor of 1.6.
Appropriate activation enthalpies for the deep Earth may induce an order of magnitude viscosity variation (q in equation 8) greater than 3. For example, assuming a CMB temperature of 3700 K, a low enthalpy of 500 kJ/mol (Yamazaki and Karato, 2001) results in q ∼ 4, and a higher estimate of 792 kJ/mol (Ammann et al., 2009) leads to q ∼ 7. When q 3, instabilities are initially confined to the boundary layer (Christensen, 1984). Further increasing q reduces the onset time to instability, decreases the scale of flow, and intensifies the time dependence (Olson et al., 1987). This eventually leads to a broad-scale well-mixed upwelling with low viscosity (e.g., Christensen, 1984;Thompson and Tackley, 1998). Our models do not generate this type of instability because of the reduced viscosity contrast (q = 3) and the imposed long-wavelength thermal perturbation. We would expect such a plume to modify the quantity of residual material at the CMB, although this regime may not be applicable to ULVZs that are only around 2 orders of magnitude weaker than the ambient material (Hier-Majumder and Revenaugh, 2010).
To associate aggregate chemical density anomalies to P-wave velocity (V P ) and S-wave velocity (V S ) of the ULVZs we use a Voigt-Reuss-Hill (VRH) mixing model (see Watt et al., 1976). We estimate ULVZ properties as a mixture of a silicate component and a dense oxide, and density anomalies are generated by varying the fraction of oxide in the solution. In our first mixing model we assume a disequilibrium assemblage to avoid unnecessary speculation on lowermost mantle mineralogy, and focus on the combination of properties rather than specific stable coexisting phases. The composition and properties of the silicate component are not well constrained in a chemically-distinct ULVZ, so we use PREM as the model for its density and velocity. The dense oxide is assumed to be (Mg 0.16 Fe 0.84 )O magnesiowüstite (Mw), as this is the only composition for which extremely low sound velocities have been measured at appropriate pressures (Wicks et al., 2010). We refer to this model as "Mw + PREM". Calculations are made at 121 GPa, the highest pressure reached in the experimental study, and our model is extrapolated to a high (4700 K) and low (2700 K) temperature estimate of the CMB. The temperature derivatives of velocity used are ∂V P /∂T = −4.64 × 10 −4 (km/s) K −1 and ∂V S /∂T = −3.85 × 10 −4 (km/s) K −1 , measured on MgO using a multi-anvil apparatus up to 24 GPa and 1650 K (Kono et al., 2010). The measured thermal expansion of FeO at 93 GPa, α = 0.54 × 10 −5 K −1 (Seagle et al., 2008), is used to estimate the temperature dependent density of (Mg 0.16 Fe 0.84 )O. Under these assumptions, our VRH mixing model produces a range of V P and V S for a given chemical density anomaly (Fig. 7A, B).
We also construct a second mixing model to consider the effects of Fepartitioning between magnesium silicate perovskite (Pv) and Mw ("Mw + Pv"). Estimates for the Fe-partitioning between these phases approaching CMB conditions varies between K Pv/Mw D ≈ 0.42 for pyrolitic composition (Murakami et al., 2005) and K Pv/Mw D ≈ 0.07 for an Al-free system (e.g., Sakai et al., 2009;Auzende et al., 2008). Trends indicate that K Pv/Mw D is likely lower at higher pressures and Fe-content, but it is not possible to reliably extrapolate from present datasets. This mixing model necessitates a small K Pv/Mw D to produce an equilibrium phase assemblage with a bulk density comparable to PREM and thus density anomalies (equivalently B) that are contained within the range of this study. Assuming (Mg 0.16 Fe 0.84 )O as the oxide phase and with K Pv/Mw D = 0.07 predicts a coexisting Pv of composition (Mg 0.72 Fe 0.28 )SiO 3 . Such high Fe -contents may be unstable in both Mw (Dubrovinsky et al., 2000) and Pv (Mao et al., 2005), but other experimental results indicate otherwise (Tateno et al., 2007;Lin et al., 2003). Possible coexistence of post-perovskite in the presence of Al and Fe adds further complication (Andrault et al., 2010), so we extrapolate assuming a high CMB temperature estimate (T CM B = 4700 K) to promote Pv as the stable silicate phase. Although the presence of Al throughout the lower mantle may alter the partitioning behavior (e.g., Irifune, 1994) and select thermoelastic properties (e.g., Jackson et al., 2005), we focus on an Al-free system to simplify the results. The composition and properties of Pv are determined from a finite strain model (Li and Zhang, 2005) where effects of Fe are neglected except for a correction to the density. In exploratory models we perturbed the bulk and shear moduli due to Fe-content (Kiefer et al., 2002) but discovered this had marginal impact. We select an adiabat foot temperature for Pv that reaches 4700 K at 121 GPa and extrapolate the dense oxide to high temperature using the same procedure as previously described.
The chemical density anomalies of the mixing models can be directly associated with the buoyancy numbers used in the geodynamic calculations (Fig. 7C). Assuming a high (1500 K) and low (500 K) estimate for the temperature drop at the CMB, we can therefore predict a relationship between the wavespeed reduction and relief (thickness) of ULVZs. Our geodynamic models assume ∆T = 1500 K (see Table 1), so we introduce an inconsistency by mapping the data using a reduced ∆T = 500 K, but the Rayleigh number is small (a factor 2-3 difference) and does not alter the dynamics significantly. Exploring the implications of a high and low temperature drop at the CMB adds to the predictive capabilities of our synthesized model. For visual clarity we now plot only the Hill bound, and compare to seismic data (Figs. 8 and 9). For example, in Fig. 8B, the "Mw + Pv" model at 4700 K (∆T = 500 K) with relief of 95.7 km produces a V P velocity reduction of ∼5% with ∼4 vol. % of Mw. In the same model, ∼9% Mw produces a V P reduction of ∼9%, resulting in a relief of 27.5 km.
A negative correlation between ULVZ thickness and wavespeed reduction is expected given the inherent tradeoff in modeling seismic data. This is more evident with V P (Fig. 8A) than V S (Fig. 9A). However, this difference may be due to sampling biases. To discern V P at the base of the mantle the original seismic studies utilize long-period data which is more sensitive to long wavelength structure. In contrast, shorter period reflected phases are commonly used to determine V S structure, and hence thin ULVZs are often detected. For both wavespeeds, a cluster of datapoints populate a broad range of wavespeed reductions for the thinnest ULVZs, implying substantial local variations.
The modeled V P and V S can reproduce the seismic observations for a reasonable range of expected temperatures and temperature drops in the CMB region. 2-D studies are denoted by filled symbols in Figs. 8A and 9A. As previously noted, the seismic data exhibits a broad scatter. Localized variations in ULVZ structure and material properties, caused by temperature or chemistry, may explain this spread. However, it is noteworthy that our synthesized model, which combines empirically determined wavespeeds with geodynamic simulations, is capable of explaining the current observations. The PREM mixing model ("Mw + PREM") intrinsically overestimates the density and wavespeeds of the silicate component in the cases we explore here, because we do not explicitly model the silicate at high temperatures (Figs. 8 and 9). This model therefore underestimates wavespeed reductions and can be considered a lower bound. The Fe-partitioning model ("Mw + Pv") provides insight into the importance of both Fe-content and temperature in modulating the physical properties and seismic characteristics of ULVZs. Since both Mw and Pv have high Fe-contents the compositional density of the assemblage is high relative to PREM and is therefore likely to form a flat layer. Higher temperatures, however, reduces total density whilst further reducing wavespeeds. At T CM B = 4700 K the density of the silicate is sufficiently reduced that a large volume fraction of Mw (also of reduced density) is incorporated into the assemblage to obtain PREM-like densities. This allows large wavespeed reductions for relatively small density anomalies.
There are multiple avenues for future work. Limited constraints on the chemistry and properties of the ULVZ phases under appropriate conditions introduce uncertainties into our VRH mixing model to predict expected V P , V S , and density anomalies. A better constraint on the temperature and phase relations in D , coupled with sound velocity measurements of these phases at CMB conditions, would help to constrain the composition. Our geodynamic simulations assume a priori that there is an availability of FeOenriched material at the base of the mantle to form ULVZs. It would therefore be useful to unravel the mechanism for creating an FeO-enriched layer to constrain its likely volume and iron content. Furthermore, our results are sensitive to the initial thickness (proxy for volume) of this layer.

Summary and conclusions
We explore the geodynamic implications of a solid iron-enriched (Mg,Fe)O layer for the origin of ULVZs. Our numerical simulations produce ULVZs with a sharp and concave-down seismic top, appropriate horizontal and vertical lengthscales and an intrinsic density increase. Models with similar features can explain the scattering origin of PKP precursors (Wen and Helmberger , 1998a). The morphology and aspect ratio (relief/half-width) of the structures are strongly dependent on the buoyancy number B (or equivalently, the chemical density anomaly) and the relief (thickness) is proportional to the initial volume of ULVZ material. Complexities in the thermal and chemical structure can promote regional variation in ULVZ character. Wicks et al. (2010) have recently demonstrated that addition of small amounts of iron-enriched (Mg,Fe)O mixed with silicates may explain the observed reductions in wavespeeds for ULVZs. We combine our geodynamic results with a Voigt-Reuss-Hill mixing model and compare with seismic data. For reasonable estimates of the CMB temperature and temperature drop, our synthesized model satisfies ULVZ wavespeed reductions and thicknesses inferred seismically, thus lending support to a solid-state origin for many ULVZs.  : (A) Seismically inferred ULVZ thicknesses and S-wave velocity reductions found globally. Following the nomenclature of Thorne and Garnero (2004): Precursors to core reflected phases Garnero and Vidale, 1999;Reasoner and Revenaugh, 2000;Havens and Revenaugh, 2001;Avants et al., 2006;Hutko et al., 2009;, travel time and waveform anomalies Helmberger et al., 1998;Wen and Helmberger , 1998b;Helmberger et al., 2000;Simmons and Grand , 2002;Rondenay and Fischer , 2003;Thorne and Garnero, 2004;Zhang et al., 2009). Filled symbols denote 2-D studies. (B, C, and D) Models using the Hill bound, for various d ch , temperature at the CMB (T CM B ), and temperature drop (∆T). Solid lines are the "Mw + PREM" model and dotted lines are the "Mw + Pv" model.