Investigating the seismic structure and visibility of dynamic plume models with seismic array methods

region, indicating that D (cid:2)(cid:2) reﬂections may not always be due to a phase transformation. We suggest that slowness-backazimuth analysis may be a useful tool to locate mantle plumes in real array data sets. However, it is important to analyse the data at different dominant periods since, depending on the width of the plume, there is probably an optimum frequency band at which the plume is most visible. Our results also show the importance of studying the incoming energy in all directions, so that any apparently out-of-plane arrivals can be correctly interpreted.


I N T RO D U C T I O N
On a large scale, convection in the mantle leads to the production of new oceanic floor at divergent plate boundaries as well as subduction of the lithosphere at convergent plate boundaries. Volcanoes situated along plate boundaries can therefore be explained by plate tectonic activities. However, surface observations such as intraplate volcanism and hotspots require a different explanation. Hotspots, such as Hawaii and Iceland, are large topographic swells associated with elevated mantle temperatures, a distinct geochemistry, and unusually large amounts of volcanism (e.g. Sleep 1990;Putirka 2008). Morgan (1971) suggested that hot, thin plumes exist in the Earth's lower mantle as a secondary mode of convection and feed the hotspots with upwelling material coming from the core-mantle boundary (CMB). Numerical modelling (e.g. Loper & Stacey 1983;Farnetani & Richards 1994;Lin & van Keken 2006) and laboratory experiments (e.g. Griffiths & Campbell 1990;Lithgow-Bertelloni et al. 2001;Schaeffer & Manga 2001;Davaille et al. 2002) have investigated the feasibility of this hypothesis.
Early studies proposed that the D layer is the likely origin for mantle plumes (e.g. Olson et al. 1987;Davies 1988). An alternative view is that plumes can also originate from a thermal boundary at the 660 km discontinuity (Schubert et al. 2001), or for some plumes perhaps no thermal boundary is required-those are sometimes called 'splash-plumes' (Davies & Bunge 2006), similar to return flow. Another hypothesis is that plumes do not exist at all, because hotspots may be caused by lithospheric melting (e.g. Anderson 2000;Foulger 2007). A way to reconcile these different hypotheses is that perhaps hotspots have a mixed origin, with only some of them being caused by deep mantle plumes (Courtillot et al. 2003).
The extent to which plumes should be purely thermal or rather thermochemical upwellings is also debated (e.g. Anderson 1975;Farnetani & Samuel 2005;Davies & Bunge 2006). This question is important because either possibility requires a different style of dynamics inside the Earth. If plumes arise from a purely thermal boundary layer, and that boundary layer is situated just above the CMB, the resultant excess temperatures inside the plume will be very large throughout the mantle (Farnetani 1997). These high temperatures may be incompatible with mantle temperatures inferred from ocean island basalts and surface heat flow (e.g. Sleep 1990;Falloon et al. 2007). A dense thermochemical layer at the base of D , on the other hand, may lead to plausible excess temperatures within the plume (Farnetani 1997). Recent numerical modelling has shown that a thermochemical plume beneath Hawaii is more consistent with images from seismic tomography than a purely thermal plume (Ballmer et al. 2013) due to the dynamic effects of chemical heterogeneity. Additionally, simulations of thermochemical plumes sourced at the base of the mantle from 'primordial' piles (Fe and Si enriched) show entrainment of dense, undegassed material in quantities that are compatible with the lowest values of the 4 He/ 3 He ratio observed in Ocean Island Basalts (OIBs) (Deschamps et al. 2011). Finally, simulations have also shown that, in addition to undegassed material, plumes may re-entrain upwards small fractions of Mid-Ocean Ridge Basalts (MORBs) injected in the lowermost mantle by slabs (Li et al. 2014), thus explaining the large dispersion in the helium isotopic ratio seen in OIBs.
Observationally, seismology can access the deep interior structures of the Earth in situ and hence, could provide important insights on the existence, nature and structure of mantle plumes. We expect plumes to be seismically slow (e.g. Goes et al. 2004;Styles et al. 2011a) because they carry hot, buoyant material, and seismic waves generally travel more slowly in such material. Measurement and inversion of seismic travel-times (i.e. tomography) is commonly used in regions where mantle plumes are expected, and seismically slow structures have been observed beneath many hotspots (e.g. Montelli et al. 2004;Zhao 2004;Montelli et al. 2006;Wolfe et al. 1997Wolfe et al. , 2009Rickers et al. 2013;French & Romanowicz 2015). While these are often interpreted as mantle plumes, various data and numerical limitations mean there are still no strong constraints on the width and depth extent of these structures (de Wit et al. 2012), which makes it difficult to infer the dynamic mechanism behind their origin.
It is especially challenging to image narrow plume stems in the lower mantle (Nataf 2000;Hwang et al. 2011;Rickers et al. 2012) due to an effect known as 'wave front healing' (Nolet & Dahlen 2000), a diffraction phenomenon which can lead to the loss of the signal from objects that are small relative to the path length of the waves (Wielandt 1987;Malcolm & Trampert 2011). Destructive interference between waves which travelled directly through a plume and waves that diffracted around it will diminish the plume signal in the first-arriving wavelet at a receiver, along the path of the direct wave. Since it is the first arrival which is typically used in seismic tomography, the plume tail becomes 'invisible'. This leads to a debate whether the plume (tail) does not exist in the first place, or whether it exists but is simply unresolvable. Additionally, plumelike structures in the upper mantle may be artificially elongated in tomographic images by vertical smearing, meaning that we may wrongly infer the presence of a deep mantle plume from a shallow structure (e.g. Maguire et al. 2018). Tomographic inversions based on fitting the full seismic waveform, rather than the traveltimes of the first-arriving phases, could potentially ameliorate these effect and allow narrow plume tails to be resolved in the deep mantle (Rickers et al. 2012(Rickers et al. , 2013. However, full-waveform inversions are time-consuming and computationally expensive (e.g. Cobden et al. 2015b), especially for the deep mantle.
Few seismic studies have presented evidence for mantle plumes with methods other than tomography. A notable exception is the study by Nataf & VanDecar (1993), who showed that systematic deviations in P-wave travel times are compatible with a plume-like structure at 700 km depth beneath the Bowie hotspot. Additionally, localised thinning of the transition zone beneath hotspots, as inferred from deflections of the seismic discontinuities at 410 and 660 km, may provide indirect evidence of thermal plumes (e.g. Shen et al. 1998), although the behaviour of seismic discontinuities in the presence of plumes can be variable and complex (e.g. Saki et al. 2015;Jenkins et al. 2016). More recently, modelling of S waveforms has revealed cylindrical ultra-low-velocity zones (ULVZs) and low velocity regions at the base of the mantle beneath Iceland and Hawaii (Cottaar & Romanowicz 2012;He et al. 2015;Yuan & Romanowicz 2017), but it is unclear if or how these deep structures may be connected with the surface hotspots.
In this study we consider the alternative possibility of using array seismology to detect mantle plumes. Array seismology is a powerful technique for the interpretation of seismic recordings (Rost & Thomas 2002, 2009. Data from several receivers placed closely together are stacked so that coherent signals will be enhanced. This way it is possible to obtain information on the directionality of the incoming energy and to detect arrivals with weaker amplitudes such as scattered signals. Scattering of seismic waves may be observed when the impedance contrast of an object with the background mantle is sufficiently large. In a numerical simulation, Dalkolmo & Friederich (2000) showed that thermal plume-like structures generate weak scattering which would be hard to detect on an individual seismogram. However, they suggested that stacking techniques could enhance the visibility of the scattered energy. Here we will analyse full seismic waveforms from simulations of wave propagation through plume models to look for scattered energy, using the stacking techniques of array seismology. This is a novel approach for investigating the size and shape of narrow mantle plumes, and may provide additional constraints on the location of plumes than can be obtained from tomography alone.
Several earlier studies have investigated scattering from synthetic mantle plumes using different methods (Ji & Nataf 1998;Tilmann et al. 1998;Dalkolmo & Friederich 2000;Capdeville et al. 2002). In these studies, the plumes were represented as vertical cylinders with sharp boundaries (or edges that decay over a narrow Gaussian window). The impedance contrasts were typically estimated from an assumed temperature contrast and seismic sensitivity to temperature perturbations. Such plume models are potentially unrealistic in that the thermal anomaly of a plume may decay over a broad distance (i.e. several hundred kilometres); the shape of the plume is more complex than a vertical cylinder; and the seismic properties of the plume are also influenced by chemical heterogeneity. In this study we use numerical modelling to define thermochemical plumes whose shape and structure are dynamically plausible. We convert these thermochemical plumes into seismically equivalent structures via thermodynamic modelling (Connolly 1990;Connolly 2005) of a recent compilation of mineral elastic parameters (Stixrude & Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Lithgow-Bertelloni 2011). We then use a spectral-element modelling code which solves the seismic wave equation (Fichtner & Igel 2008) through the plume models, and generates synthetic seismograms for receivers at the surface. By solving the seismic wave equation numerically, scattering effects are implicitly included and no approximations of scattering behaviour (as was typically needed in earlier studies) are required.
The receivers are divided into several arrays and we process the composite seismograms of each array using methods of array seismology. This is a feasibility study intended to establish how plume-like structures might be observed with array methods.

B A C KG RO U N D : A R R AY M E T H O D S
Array seismic recordings produce homogeneous data sets that can be used to study the fine-scale structure of the Earth's interior (e.g. Rost & Thomas 2002,2009). One of the key advantages of array data sets, which we will exploit in this study, is that directionality of the incoming signal can be estimated.
Array analysis is based on the assumption (approximation) that the incoming energy from a seismic event behaves like a plane wave. A plane wave arriving at an array of seismic stations will have a different arrival time at each station (delay time), depending on the slowness (1/wave speed) of the wave and the geometry of the array (Schweitzer et al. 2002). Only for stations aligned parallel to the incoming wave front can no time-shift be recorded. There are many techniques in array seismology; pertinent for this study are Velocity Spectral Analysis (in short: Vespa process (Davies et al. 1971)) and slowness-backazimuth analysis (King et al. 1976;Rost & Thomas 2002). The former is useful for investigating either the slowness or backazimuth of the incoming energy as a function of time. The latter is useful for studying slowness and backazimuth simultaneously in a given time window.
In the Vespa process, the delay time corresponding to a chosen (3-D) slowness vector is subtracted from each seismic trace, before the shifted traces are summed together (i.e. stacked). This is called the 'delay and sum' or beamforming method, and the result is that uncorrelated noise and phases with different slowness values than the chosen one are suppressed, while phases with the prescribed slowness and backazimuth are amplified. When this procedure is done for a range of slowness vectors, the stacked traces can be displayed over time as a function of the horizontal slowness for a fixed backazimuth, or vice versa. The resulting diagram is called a vespagram.
We use 4th root vespagrams for improved slowness resolution (Muirhead & Datt 1976;Rost & Thomas 2002). This involves taking the 4th root of each seismic trace before stacking them. The beam trace, i.e. the resultant time-series after stacking, is then raised to the 4th power. 4th root vespagrams boost the amplitude of small coherent signals relative to larger ones, hence enhancing the visibility of the former. However, because the waveforms are distorted then the vespagrams cannot be used for interpretations regarding amplitude or polarity. An alternative to the 4th root vespagram is to make a phase-weighted stack (Schimmel & Paulssen 1997). For moderate-sized arrays the differences between phase-weighted and 4th root vespagrams are negligible (Rost & Thomas 2009). In this study we use vespagrams of slowness versus time. This does not directly tell us about scattered energy but helps us to identify which seismic phases we are looking at in a given time window.
For studying scattering, it is important to know what is theoretically the backazimuth for energy travelling along the great circle path between source and receiver. Usually this is calculated from a reference station located in the geographical mid-point of the array. We can then study the energy arriving at the array as a function of slowness and backazimuth within a given time window. This slowness-backazimuth analysis is described in detail in Rost & Thomas (2002), and is similar to frequency-wavenumber (f-k) analysis, but the visualisation of results is different (Rost & Thomas 2009). When a wave travels along the great-circle path from source to receiver, it will have the same backazimuth as the theoretical backazimuth. However, if a wave encounters a heterogeneity, the wave fronts are distorted and seismic energy can be deflected. This means that energy arrives at the receivers with a different slowness and backazimuth than the theoretical one. Those arrivals are called 'out of plane' signals (OOPs). For interpretation we typically assume that an OOP wave has been scattered once between source and receiver (waves that have undergone multiple scattering would probably have lower amplitudes). The observed slowness of the signal depends on the depth at which it was scattered and is constrained by the distances between source, receiver and obstacle. In theory, OOPs can therefore be back-traced and the location of the scattering object can be determined. This method has been applied with success to locate subducted slabs (e.g. Kaneshima & Helffrich 1999;Rost et al. 2008;Weber et al. 2015;Schumacher & Thomas 2016) and large low-velocity regions (Schumacher et al. 2018), and we are curious if it could also be used to locate mantle plumes.

Dynamic plume models
Our plume models are generated from numerical simulations of mantle convection in spherical coordinates, using the code StagYY (Tackley 2008). This code solves the equations of conservation of mass, momentum and energy for a compressible fluid with infinite Prandtl number. Thermochemical models further require solving an equation for the conservation of composition. Composition is modelled with a collection of tracers of two types representing the dense and regular material, respectively, with dense tracers being initially distributed in a basal layer. At each timestep, the compositional field is calculated from the concentration of dense and regular tracers in each cell. Note that the exact chemical nature of the compositional field is not prescribed a priori-it is described only in terms of a density contrast between dense and regular material-and can thus be decided a posteriori. We restrict our calculations to a sector with dimensions π 4 × π 8 (i.e. 45 • × 22.5 • ) in latitude and longitude. In the vertical direction our models span the whole mantle down to the CMB. The models are computed on a grid with 128 × 64 points horizontally, and 128 points vertically, with grid refinement at the top and bottom, allowing a good description of the thermal boundary layers. Grid refinement is also prescribed around the 660 km boundary, where plumes are impeded and thinned due to the combined effects of the ringwoodite to perovskite phase transition and of the viscosity change (here fixed to a viscosity ratio of 30 between the lower and upper mantles (e.g. Mitrovica & Forte 2004). The boundaries at the surface and the bottom of the model are free slip (i.e. the horizontal stress is zero) and isothermal. The sides of the domain are also free-slip and reflecting, which means that the solution calculated within the domain may be extended outside this domain by performing a planar symmetry with respect to each side. The system is heated from below and internally, with the latter being equivalent to a surface heat flow of ∼10 mW m -2 . Further details of the convection modelling can be found in the Appendix and Deschamps et al. (2018). Note that in simulations the output temperatures ( Fig. A1) are in non-dimensional form and uncorrected for the adiabatic effects. For our purpose, we re-dimensionalized temperature and added an adiabatic effect following: here T S = 2500 K is the superadiabatic temperature jump; a(z) is the adiabatic correction at depth z; T(x, y, z) is the non-dimensional temperature at a given location, andT top is its surface value, which we fixed to 0.12, a value that is equivalent to a dimensional surface temperature of T surf = 300 K (see the Appendix).
The flow structure, including the size and shape of plumes, depends very much on the choice of input parameters for the numerical simulation, most particularly on the chemical density contrast ( ρ C ) between dense material and regular mantle (e.g. Davaille 1999;McNamara & Zhong 2004;). Here, we consider two different models for two chemical density contrasts (measured with the buoyancy ratio, B, which describes the ratio of the chemical density contrast ( ρ C ) to the thermal density contrast, see the Appendix, eq. A1): B = 0.23 (corresponding to a density contrast ρ C = 140 kg m -3 ) and B = 0.18 ( ρ C = 110 kg m -3 ), which we shall refer to as TC1 and TC2, respectively. These buoyancy ratios were chosen to ensure that the reservoirs of dense material at the base of the plumes remained stable over time, and correspond to an enrichment in iron of 2-4 per cent, which is consistent with the ranges inferred from seismic studies of LLSVPs (e.g. Trampert et al. 2004;Mosca et al. 2012). Except for the chemical density contrast, all other properties of TC1 and TC2 are similar (see Table A1). The thermochemical plumes obtained in these two models are substantially different (Figs 1 and A2). In both cases the plumes are generated by thermal instability at the top of a chemically dense pile. However, despite the chemical origin for initiating the plumes, they are predominantly thermal anomalies, with only a narrow trail of chemically heterogeneous material being entrained along the central axis of the plume. When plotted in terms of absolute temperature ( Fig. 1), TC2 has a wider plume stem than TC1 but rises from a more compact base. Globally, model TC2 is also hotter than TC1. However, when the plume boundary is defined in terms of a thermal anomaly (Table 1, Fig. A2) the thermal anomaly of the plume itself in TC2 is not significantly hotter than in TC1, and the widths are comparable. More importantly, the plume in TC2 entrains a larger amount of dense (i.e. chemically distinct) material. The vertical advection of material is impeded by the viscosity change at the 660 discontinuity, so that both plumes spread out below the discontinuity, and then thin significantly after rising through it. Impingement of the plume head at the surface causes the plume to spread outwards and then downwards, leading to perturbations in the temperature and chemistry at distances far (>1000 km) from the plume axis. Also noteworthy is that the plumes generated in both TC1 and TC2 are affected by slight tilting due to subhorizontal flow. Downwellings of material are located on the sides of the model, a distribution that is constrained by the geometry of the system.
Defining a plume diameter is somewhat subjective, since the thermal anomaly decays over a large distance (see Fig. A2), and may not coincide with the resultant seismic anomaly. Using Fig. A2 we estimated the width of the plume at four representative depths, summarised in Table 1. Our criterion for the plume boundary is the location where the temperature anomaly has decayed to 50 per cent of its maximum. Alternative criteria, such as the location of the maximum temperature gradient, would give different dimensions.
Prior to conversion to seismic structure, the plume models are recalibrated in the vertical direction so that the grid points are uniformly spaced every 15 km. This results in 193 gridpoints and is done so that the models are appropriately formatted for the seismic modelling code.
In numerical simulations, the compositional field C varies between 0 for cells filled with regular material (tracers) only and 1 for cells with dense material only. The nature of regular and dense material is not prescribed a priori and can be fixed by the user. We assume that the regular composition (C = 0, C 0 ) is a pyrolite (Sun 1982) with the following bulk composition (molar percentages): 49.13% MgO, 38.61% SiO 2 , 6.24% FeO, 3.25% CaO, 2.77% Al 2 O 3 and that dense material (composition C = 1, C 1 ) corresponds to a 4% increase in the molar percentage of FeO compared to pyrolite, and has the composition: 47.1% MgO, 37.0% SiO 2 , 10.2% FeO, 3.1% CaO and 2.6% Al 2 O 3 . This choice for dense material is motivated by the observation that LLSVPs might be enriched in iron (Trampert et al. 2004;Mosca et al. 2012).

Conversion from thermochemical to seismic structure
We convert our 3-D thermochemical plume models into 3-D seismic models using the thermodynamic modelling code Perple X (Connolly 1990(Connolly , 2005. For a given pressure P, temperature T and composition C, Perple X uses a set of mineral thermodynamic and elastic parameters, together with an equation of state, to compute the equilibrium mineral assemblage and associated seismic properties. Since our plume models are defined as a function of depth rather than pressure, we must convert all the depths to pressures, and we use PREM (Dziewonski & Anderson 1981) for this calibration. Within Perple X, we use the equation of state of Stixrude and Lithgow-Bertelloni (2005), and the mineral parameters of Stixrude & Lithgow-Bertelloni (2011). The average seismic properties (density, P-wave speed and S-wave speed) are found by computing the Voigt-Reuss-Hill average for the mineral assemblage at each P, T and C.
This procedure ensures consistency between the seismic and thermochemical properties of the plumes, especially since the mineral elastic properties and phase boundaries are defined within the same thermodynamic framework. Our method is similar to other recent studies, in which dynamic plume models were converted to seismic structures (e.g. Hwang et al. 2011;Styles et al. 2011a;Maguire et al. 2016Maguire et al. , 2018. Those earlier studies were primarily investigating the influence of plumes on seismic traveltimes and waveforms. The main difference between our procedure and these other studies is that here the plume dynamics are computed with full thermochemical convection, rather than purely thermal convection. While our plumes do not entrain much chemical heterogeneity ( Fig. 1), allowing for such variability is more consistent because the details of the flow and thermal fields are strongly affected by the presence of chemical heterogeneities. In particular, it changes the temperature anomalies required to generate the plumes (e.g. Farnetani 1997) and potentially some details of the thermal structure as they ascend.
The seismic wave speeds of our plume models have not been corrected for anelasticity. Anelasticity causes a reduction in the wave speeds at high temperatures. In the lower mantle, this is likely of negligible consequence (Brodholt et al. 2007), but in the upper mantle it means that we may underestimate the wave speed anomaly Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 of the hottest parts of the plume relative to the surroundings. As this is a feasibility study and we are not trying to quantify the exact temperatures at which plumes can be detected, this effect is unlikely to influence our conclusions, and in fact means the results are more conservative. The seismic properties of the plume models (wave speeds and density) are shown in Fig. 2.

Reference model
In seismic tomography, 3-D seismic structures are usually shown as deviations from a 1-D reference model rather than absolute wave speeds, in part because the magnitude of lateral variations are typically only a few per cent. Commonly used reference models are PREM (Dziewonski & Anderson 1981) and AK135 (Kennett et al. Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020  Fig. A1, with the criteria that the plume boundary temperature is given by where <T> is the average temperature at a given depth, T max is the maximum temperature at that depth, and c is a constant. Here we assume c = 0.5.

Depth (km)
Width of thermal anomaly (km) Amplitude of thermal anomaly (K)  TC1  TC2  TC1  TC2 500 1995) which are derived from global seismic data set and can be thought of as the average 1-D seismic structure of the Earth. For our interpretation purposes, it is also useful to have a 1-D 'no plume' reference model, so that we can be certain which signals in our data are due to the plume, and not some other factor. Additionally, since the wave speed and density variations associated with the plume are subtle (Fig. 2), it is helpful to enhance their visibility by showing the seismic anomalies relative to a reference model. We tested five types of 1-D reference model: (i) the seismic references AK135 and PREM; (ii) the horizontally-averaged seismic structure of model TC1; (iii) the horizontally-averaged thermochemical structure of TC1, converted into seismic structure; (iv) a vertical profile at the edge of model TC1, far from the plume centre and (v) the seismic structure corresponding to pyrolite, i.e. composition C 0 , with an adiabatic temperature gradient, and range of potential temperatures.
In each of models (i)-(iv), we found that the seismic difference between the plume models and the reference is dominated by differences in discontinuity structure between 410 and 660 km. These differences are amplified by the recalibration to 15 km spacing between depth points. The differences span the entire model and are so strong that they overprint the seismic anomaly of the plume (see Figs A3 and A4). This makes it impossible to distinguish between seismic effects of the plume versus the peculiarities of the reference model. Seismic reference models such as PREM and AK135 are the seismically sampled 1-D average of the Earth's interior, and have been shown to not always correspond to a physically meaningful structure, especially in the vicinity of the transition zone (e.g. Cobden et al. 2008;Styles et al. 2011b). Likewise, models (ii) and (iii) both involve taking a weighted 1-D average of the 3-D structure, which can lead to non-physical discontinuity depths. Meanwhile, the limitation of reference (iv) is that all points in the model are participating in the convection and, due to the small size of the box, are influenced by the plume, so there is no 'plume-free' part of the model. In all these cases, we cannot reasonably ascribe the differences in seismic structure between the reference and plume models as an effect of the plume as opposed to the choice of non-physical reference model. However, we found that for model (v), the adiabatic pyrolite, with a potential surface temperature of 1473 K, the depths of the 410 and 660 discontinuities are mostly similar to those in the plume models, except near the central axis of the plume, and hence most of the differences in seismic properties are related to the plume conduit (see Fig. A4). We therefore use this as our 'no plume' reference model.
The reference model and its seismic properties are shown in Figs 1 and 2, respectively. The differences between the reference model and the two plume models are shown in Fig. 3. We emphasise that the purpose of Fig. 3 is to enhance the visibility of the plumes. The magnitudes of the seismic anomalies are dependent on the choice of reference model and should not be interpreted as representative of what may occur in the real Earth or in seismic tomography models. In Fig. 4 we compare 1-D vertical profiles through the plume and the reference models, with profiles taken through the centre of the plume and far from the centre.

Waveform modelling
To simulate waves travelling through the different plume models we use the code SES3D (Fichtner & Igel 2008;Gokhberg & Fichtner 2016). This uses a spectral-element solver for the calculation of the forward wavefields (for details see Kopriva 2009). The integral form of the elastic wave equation in a heterogeneous medium is solved on a hexahedral grid in spherical coordinates. It is convenient for the purpose of this study because it is designed for regionaland continental-scale simulations and hence suitable for capturing the dimensions and internal details of mantle plumes. Unwanted reflections from the model boundaries are mostly suppressed due to absorption by a perfectly matched layer (PML).
The whole model domain spans 45 • in latitude, 22.5 • in longitude, and 2615 km in the vertical direction. For computational simplicity, SES3D does not contain a CMB and since we are only interested in the effects of the plume conduit such a set-up is reasonable.
We run two sets of simulations. For the first set, the computational grid is divided in 96 × 48 × 50 elements in latitudinal, longitudinal and vertical directions. Consequently, the maximum size of each element is roughly 50 km in each direction. Within each element, the wave equation is discretised on 125 grid points. One wavelength should be at least as long as the size of 1.5 elements to ensure accurate solutions, which in theory leads to a minimum wavelength of λ = 75 km. We used a Heaviside function filtered with band pass between 25 and 100 s as a source-time function, so the dominant period is 25 s. This roughly translates to S-wavelengths of 130 and 180 km in the upper and lower mantle, respectively. The corresponding P-wavelengths are approximately a factor of 2 larger.
In the second set of simulations, the computational grid is divided into 162 × 84 × 82 elements in latitudinal, longitudinal and vertical directions. This yields 30 km as the maximum element size and so allows shorter wavelengths down to λ = 45 km for accurate computations. The source-time function in this simulation has a dominant period of T = 15 s.
The sampling rate is chosen to be 0.09 and 0.03 s for the lowerand higher-frequency simulations respectively. For both simulations the recordings are roughly 1200s long. The source has an explosiontype moment tensor as focal mechanism and is placed in the middle of the short side of the model (see Fig. 5) at a depth of 500 km. Therefore, interference of direct waves and depth phases can be excluded. The position of the plume is roughly 20 • epicentral distance from the source. On the surface, receivers are evenly distributed with an interstation spacing of 2 • , as can be seen in Fig. 5, giving a total of 180 seismic stations.

Array analysis
We divide the 180 receivers into 15 arrays of 9 to 16 stations, as shown in Fig. 5. Since the stations are spaced 2 • apart, the apertures of the arrays are between ∼4 • and 6 • . Some of the stations are used in more than one array. Because the arrays are rectangular, there is usually no station at the geographic centre of the array to take as Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 the reference station. We therefore take the station which is second from the left and the bottom in each array (relative to the plot in Fig. 5) as the reference and calculate the theoretical backazimuth from this station.
We compute vespagrams for arrays 1-8. The epicentral distances range from 32 • to 39 • , and this limits the maximum depths of the direct P and S wave to ∼1020 and 1000 km, respectively. The four arrays on the left side of the model (as plotted in Fig. 5), 9-12, are used to investigate possible backscattering, i.e. energy that is reflected in a very steep angle and travels back into the direction of the source. Arrays 13-15 are used to look for OOPs directly above the plume, i.e. if energy gets deflected straight upwards by the plume. In each array, we computed sliding-window slownessbackazimuth diagrams for the R, T and Z components. The time windows are 70 s long and the window is shifted in intervals of 35 s (which leads to an overlap of time windows).

R E S U LT S
The output of the wave propagation simulation is in the form of a three-component seismogram (radial = R-component, transverse = T-component, and vertical = Z-component) for each receiver placed on the surface. Example seismograms are shown in Fig. 6 for the simulations at 25 and 15 s. The seismograms are taken from the station indicated by a bold cross in Fig. 5, which is 33 • epicentral distance from the source. This station is positioned beyond the plume such that waves travelling here will have sampled the plume conduit. The arrival times of the main seismic phases as predicted by ray theory, calculated with TauP toolkit (Crotwell et al. 1999), are also indicated. The differences in the waveforms and traveltimes of the earlier arrivals for the plume and 1-D reference (no-plume) models are very small. The first arriving wavelet, the direct P wave, arrives at around 380 s with a bottoming depth of 890 km for 33 • distance, and shows almost no difference in the seismograms of the reference model (black), TC1 (blue) and TC2 (red). In the later arriving waveforms, differences in the amplitude and shape of the waveforms can be observed, although they are more pronounced for the 25 s simulation than the 15 s one.
For the array analysis, we always compare waveforms for the plume models with the same waveforms for the reference model. This enables us to distinguish which signals are influenced by the plume, and how. Ideally there should be no OOPs in the reference model as the structure is entirely 1-D. In practise, we occasionally see OOP energy in simulations for the reference model, which arises from imperfect absorption (i.e. partial reflection) of energy at the model boundaries. Numerically, it is impossible to suppress these reflections completely. Such OOP phases are weaker than the direct arrivals and may occur over a range of backazimuths. We therefore only interpret OOPs (and other phases) in the plume models as resulting from the plume(s) if the same phases are absent in or  deviated in slowness or backazimuth from the reference model. As there are many results from this investigation, we show only the key findings in the main text. Additional information can be found in Appendix II.

Array analysis for 25 s dominant period
The differences between the two plume models and the reference model are clearest in the Z component. In 11 different cases, we observed out-of-plane signals in the plume models, that are inplane for the reference model. The clearest examples are shown in Figs 7-9, together with the corresponding vespagrams (the OOP phase is circled in red). At shorter distances from the plume (Array 7) the OOP signals correspond with the expected arrival time and slowness of an sP phase; at intermediate distances (Arrays 2 and 3) to a pP phase, and at longer distances (Array 4) to an sS phase. In the vespagrams (which are calculated for the theoretical backazimuth), the waveform for the OOP phase does not vary significantly between TC1, TC2 and the reference model. Note also in the slownessbackazimuth diagrams, that arrival energies are indicated by a colour scale which is relative to the amplitude of the first arriving phase (i.e. the direct P wave), and the total energy is in arbitrary units. This means that a signal which appears stronger in one diagram than in another may not necessarily have the larger energy, if the scaling is different.
The OOPs in TC1 and TC2 are deviated by ∼5-10 • relative to the great circle path. Sometimes the arrival in TC2 is more deviated than TC1 (see Fig. 9). Interestingly, when the array is at a positive backazimuth with respect to the plume then the deviations are negative, while when the array is at a negative backazimuth with respect to the plume then the deviations are positive.
In the arrays with larger epicentral distance, we also observe a signal in the plume models that does not appear in the reference, but is travelling on the great circle path. It has a slowness of about Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 4 s deg -1 and arrives at ∼850-900 s (Fig. 10), corresponding to a phase which was reflected somewhere near 2500 km depth. The significance of this extra arrival is discussed in Section 5.2.

Array analysis for 15 s dominant period
Since the computational effort at 15 s is significantly larger than at 25 s, and since all observed effects at 25 s were stronger for TC2 than TC1, we only performed the 15 s simulation for the reference model and TC2. At 15 s we see more scattered energy, and arrivals which previously appeared as one large (smeared) phase at 25 s can often be identified as two separate phases. We sometimes see OOP signals in TC2 that are also present in the reference model, and so they cannot result from the plume. Most likely they are due to partial reflections at the boundaries of the model. Comparing the plume to the reference model, we also observe changes in the slowness and strength of waveforms with arrival time and slowness corresponding to sS and/or SS that are difficult to interpret.
Often, we see the 'array response function' (ARF) in the slowness-backazimuth plots, in other words side-lobes of energy around the main phase. It is unsurprising to observe this at 15 s, since we are now modelling waves with a higher dominant frequency. The side lobes by definition are deviated from the theoretical backazimuth, but we should only interpret the signals as actually being out-of-plane if the whole ARF is deviated relative to the theoretical backazimuth.
In the Appendices (Figs A5-A7) we show the equivalent slowness-backazimuth diagrams for the 15 s simulation, as those shown at 25 s (in Figs 8-10). At 15 s, the same arrivals that were deviated at 25 s are slightly deviated for the plume model relative to the reference, and with the same direction of rotation as at 25 s. However, because the arrivals are less clear, and their deviations are small (typically <5 • , i.e. within observational error), we do not consider them as a robust observation. We also sometimes see energy with a slowness of ∼2 s deg -1 which is OOP for both the plume and the reference model, but strongly deviated in the plume model compared to the reference (by 10-20 • ), and with the same sense of rotation as the other deviated phases. However, a slowness of ∼2 s deg -1 corresponds to reflections somewhere near the bottom of the model. While there may be a small seismic discontinuity close to the bottom of the reference model (see Fig. 4) it is unclear if the OOP signals in the reference model are due to this discontinuity, or are an artefact of the model setup (i.e. partial reflections, because the boundaries are not 100 per cent absorbing). For this reason we do not interpret these signals.

Apparent out-of-plane signals
The out-of-plane arrivals that we observe at 25 s have a backazimuth which is opposite to what we would expect if they were generated by scattering from the plume (see cartoon, Fig. 11). The backazimuths of these arrivals are, instead, consistent with the wave front 'bending' around the plume conduit, as illustrated in Fig. 11. Array seismology assumes that the incoming energy is arriving as a planar wave front. If a part of that wave front has sampled the plume and a part has not, the result is a partial retardation of the wave front (in the part of the wavefield that travelled through the plume). On stacking the arrivals for slowness-backazimuth analysis, with the assumption of a planar wave front, this may then give signals that appear to be deviated from the great-circle path, with a sense of rotation towards the plume rather than away from it (the latter would happen for scattering). The fact that we can observe this better at lower frequencies is essentially an interference effect that works to our advantage. At higher frequencies (15 s), waves that travel through the plume versus waves that travel next to it can be better distinguished from each other, and so the OOP arrivals are less obvious.
We do not see any out-of-plane signals for the arrays located above or behind (meaning, on the source-side of) the plume (Arrays 9-15, Fig. 5), i.e. no evidence for backscattering. However, Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 the arrays are relatively close to the plume, and the stations are widely spaced, such that the difference in backazimuth between the stations may be too large for array stacking to produce coherent signals. In that case, scattering may still be evident in the original seismograms.

Reflection at 2500 km depth
For the simulation at 25 s period, we occasionally detect a phase in the vespagrams for TC1 and TC2 which arrives at ∼850 s with a slowness of 4 s deg -1 (Fig. 10). The phase is not out of plane, and is not present in the reference model (see Fig. A8). Using the Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 TauP Toolkit and the 1-D seismic profile in TC2 at the theoretical bouncepoint for this phase (Fig. 10), we found that it matches the traveltime and slowness of an sP2500S phase, which is a topside reflection and a conversion of an sP wave into an S wave at 2500 km depth (P2500S would arrive ∼100 s too early).
Many observational studies have found seismic reflections at a similar depth, for example, the D discontinuity which has been observed ∼100-400 km above the CMB (e.g. Wysession et al. 1998). The D discontinuity is often attributed to a phase transition from bridgmanite to post-perovskite (e.g. Murakami et al. 2004;Oganov & Ono 2004;Wookey et al. 2005). At first glance this might be an obvious candidate for the reflection in our plume model. However, according to Perple X, post-perovskite is not stable at the temperatures and pressures (∼3100 K, ∼120 GPa) inside the plume core (the bouncepoint of the reflector), for both C 0 and C 1 . While the real position of the phase boundary is subject to many (experimental) uncertainties (Cobden et al. 2015a), this indicates that, at least within the framework of our modelling, post-perovskite is not responsible for the reflection at 2500 km.
At the same time, 2500 km corresponds to the top of the (thermo)chemical pile which serves as the plume source (see Fig. 10). The changes in seismic properties at this depth-i.e. a sharp decrease in both P and S-wave speed of a few per cent, but an increase in density (Fig. 10)-are also what we would expect for an increase in iron content (e.g. Deschamps et al. 2012). We suggest, therefore, that this is responsible for generating the reflection seen in our vespagrams, demonstrating that it is possible to generate seismic reflections from thermochemical piles, which may be present at the CMB in addition to post-perovskite. This is consistent with our previous work (Cobden & Thomas (2013), in which we found that the detailed structure of D reflections is sometimes-in certain locations-more easily explained by chemical heterogeneity than a phase transformation from bridgmanite to post-perovskite.

D I S C U S S I O N A N D I M P L I C AT I O N S
It is perhaps counter-intuitive that we see clearer signals from the plume at longer periods than at shorter periods. However, the outof-plane signals that we observe are a low-frequency phenomenon in which wave behaviour can be represented by planar wave fronts that curve around the plume conduit. At higher frequencies the wave behaviour becomes more complex and the coherency of the OOPs is reduced, leading to non-detections. The exact frequency Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Figure 11. Schematic diagram illustrating how the apparent OOP signals seen at 25 s may be generated, by bending of the wave fronts around the plume (left-hand panel), to give the observed slowness-backazimuth (sloaz) signals (right-hand panel). The plume is slowest at its centre hence the waves are more strongly deviated towards the centre.
at which the apparent OOP signals are clearest will depend on the width of the plume (and to an extent, the aperture of the array). In our particular example, 25 s works better than 15 s, for plume diameters of the order 400-500 km at 500-1000 km depth. In future studies, for example to test the method against real data, it would be important to filter the data at a range of frequencies in order to find the optimum frequency. This could then be used to estimate the diameter of the object causing the waves to bend.
With regard to scattering and back-scattering, which we do not observe, there are two possible explanations: either the impedance contrasts between the plume and the surroundings are insufficiently large to generate observable scattering, or the configuration of our model-in terms of source, plume and receiver locations, or the seismic frequencies-cannot resolve scattering. To test the former hypothesis, future studies with different synthetic plumes (different widths and temperature contrasts) could help to quantify the relationship between plume dimensions and optimum observational frequency, and establish whether scattering may appear for other types of plume. To test the latter, we should implement shorter interstation distances-which may help to identify less coherent scatteringand higher seismic frequencies. The arrays in this study are relatively large (up to ∼6 • aperture) with a coarse interstation spacing (2 • ). Both smaller aperture arrays and finer interstation spacings are more capable of capturing coherent signals from shorter wavelengths, but smaller arrays generate broader maxima in the waveforms, i.e. the errors in picking the slowness and backazimuth of the phases are larger. Hence the ideal configuration for detecting highfrequency scattering and with a high slowness-backazimuth resolution would be large arrays with many closely spaced stations. The same arrays would also be suited to identifying wave bending at the lower frequencies, and perhaps determining a frequency-dependent transition from wave bending to scattering. However, with current computational resources, it is difficult to run full-waveform simulations below about 10 s, over length-scales relevant for mantle plumes.
One of the limitations of the modelling in this study is the restriction to a box with dimensions 45 • × 22.5 • . This size was chosen to optimise the available computational resources, but it means that the direct P and S waves cannot sample deeper than about 1000 km. For more detailed investigation of the plume conduit in the lower mantle, it would be necessary to have a bigger box. This may also be a contributing factor for the non-observations of scattering-perhaps larger epicentral distances are required to sample the full range of scattered phases. We are currently working on extending the model domain so that the direct waves will sample the entire lower mantle, in addition to running simulations at higher frequencies and with shorter interstation distances.
This study can hence be thought of as a 'proof of concept' study, to probe what potential array methods have for detecting plumes. In a future, more quantitative study, the seismic wave speeds should be corrected for anelasticity. This would enhance the velocity contrasts of the plume in the upper mantle, and hence potentially its visibility also. Eventually, a more realistic plume model could incorporate anisotropy (alignment of grains due to flow in and around the plume) and partial melting, although computationally this would be non-trivial. It would also be interesting to investigate the effect of different viscosity ratios between upper and lower mantle-as this could change the amount of material ponding beneath 660 km-and different buoyancy ratios. Our study highlighted potentially interesting signals from the D region. Further investigation of these phases and also those interacting with the core, would necessitate using a waveform modelling code that includes the CMB. For full consistency, the dynamic modelling should then also incorporate the phase change from bridgmanite to post-perovskite, although global Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 convection simulations indicate that the influence of post-perovskite on large-scale structures is relatively limited (Li et al. 2015).
Some future refinements in the geodynamic modelling would also be beneficial for this sort of study. At present the depth of the 660 phase boundary in the geodynamic models-imposed a priori at 662 km, with a given Clapeyron slope and viscosity jumpmay not match exactly with the depths of the transition zone phase boundaries computed by Perple X. Hence, the dynamic properties of the plume may not be entirely consistent with the seismic properties in the vicinity of 660. However, the computational grid for the full waveform modelling, with elements of 30-50 km, is coarser than that of the geodynamic modelling. This is bigger than the discrepancy in the discontinuity depths between the two methods, and our results should not be affected by the inconsistency. For more detailed modelling in future, it would be desirable to align the phase boundaries in the geodynamic modelling with those in the thermodynamic modelling. This has been done previously for global-scale simulations (Nakagawa et al. 2009). Another issue is that the density defining the buoyancy of the plume may not be fully consistent with the density inferred thermodynamically with Perple X, due to slightly different ways in which chemical heterogeneity is defined. Although it would be useful to address this effect, doing so is beyond the scope of this study and we do not anticipate it to alter our results significantly.
The OOP phases that we do observe correspond in time and slowness with the depth phases sP, pP and sS. We note, however, that these phases partially overlap (in time and slowness) with the underside reflections PP and SS, respectively-see Fig. 6-and the two are hard to separate on vespagrams. Studying 1-D profiles from our plume models with TauP shows the OOPs are more likely to result from the depth phases than the underside reflections. However with a source at 500 km and the given epicentral distances of the arrays, there is also potential for triplication of both these types of phase, and when the waves do interact with the plume they may be deviated in both slowness and backazimuth, which makes it hard to distinguish one from the other in the data. The depth phases have a maximum depth of ∼750-800 km, and where they intersect the plume they are most likely traversing close to 660 and the top of the lower mantle. It is interesting to consider why we see the OOPs for these phases and not others, for example, the direct P and S waves. We suggest it may be related to the plume showing stronger perturbations in wave speed and density in the vicinity of 660 (Figs 2-4), and also a broadening of the central conduit (most visible in Figs 1, 3 & 4), i.e. the plume is both seismically stronger and wider in this region. The changes in seismic structure may be caused by (1) changes in the mineral phase relations at higher temperatures (see 1-D profiles, Fig. 4, e.g. Jenkins et al. 2016) and (2) the 660-discontinuity serving as a partial barrier to upward flow, causing plume material to spread out beneath the boundary.
Real data contains noise, and plumes may not be the only seismically significant structure in a region of observation. Both of these factors can reduce the visibility of the kind of signals which we see in our synthetic study. The effect of noise on slowness-backazimuth analysis has been investigated extensively in previous studies such as Kito et al. (2008) and Rost & Thomas (2009). These tests have shown that when the noise level exceeds 50 per cent of the main P-arrival, the analysis will break down. In real data we only usually work with data with a signal-to-noise ratio above 3. It is of course difficult to assess whether OOPs would still be visible in real data, without performing a study on real data. However, a related study on ULVZs (Cottaar & Romanowicz 2012), applied beamforming methods to S-wave data sampling the CMB region below Hawaii.
The authors observed a systematic shift in the backazimuth of the waves, consistent with refraction of the waves at the northern margin of the Pacific LLSVP, southwards into the LLSVP. While an LLSVP is much broader than a plume, the cause of the backazimuth deviations is similar to what we observe in our synthetic study, and the fact that systematic deviations of the order of 2-3 • can be resolved by the data lends support to the idea of using backazimuth analyses to locate mantle plumes. The advantage of studying slowness and backazimuth simultaneously (as in this study), is that it allows us to place a constraint on the likely depth of the object causing the deviations.
In terms of distinguishing signals from plumes versus other objects, it is sensible to begin investigations with source-receiver configurations that sample the mantle beneath known hotspots. Some examples would be earthquakes in South America travelling to the KNET array, which could sample possible plumes beneath the Canary Islands, or events from Tonga-Fiji arriving at the US-Array, which will have passed beneath Hawaii. Both of these configurations involve larger epicentral distances than were used in our modelling, which is why our ongoing synthetic studies at larger distances are important. For accessing plumes at shorter epicentral distances, large ocean experiments placing OBSs in circular deployments around hotspots would be extremely useful.

C O N C L U S I O N S
We performed a multidisciplinary numerical simulation in order to generate synthetic seismic array data for waves traversing a 45 • × 22.5 • region of the mantle with and without mantle plumes. Using array methods (slowness-backazimuth diagrams and vespagrams) we could not detect scattering from the plumes, but we did see out-of-plane arrivals consistent with an apparent bending of the wave fronts around the plume conduit. These arrivals are more obvious when the data is filtered at longer dominant periods (i.e. lower frequencies). Slowness-backazimuth analysis is potentially a novel tool for identifying mantle plumes in real data and constraining their diameter. It is relatively simple to apply and can incorporate large investigation areas. Because our study is limited in its spatial dimensions, we could not sample the plume conduit in detail in the lower mantle. Hence we suggest further modelling studies are warranted, as it is likely that further expression of mantle plumes will exist in array data, when larger model dimensions are used. Additionally, using denser arrays and modelling with shorter seismic periods may identify scattering or back-scattering that could not be detected with the present model set-up.
Our study highlights the importance of considering waves travelling at both longer and shorter periods, as there is likely an optimum wavelength for locating a plume. Furthermore, since the apparent azimuth of out-of-plane signals is different depending on the cause of those signals (e.g. wave scattering versus bending), then in real data we would need to consider waves travelling in all directions, in order to correctly interpret those signals.

A C K N O W L E D G E M E N T S
The author contributions are as follows: This paper was developed from a master's thesis by FS. FD performed the geodynamic calculations in StagYY and provided associated figures. All other aspects of numerical modelling and data analysis, including the seismic waveform modelling in SES3D, and array analysis using Seismic Handler (Stammler 1993), were done by FS. FS supplied most of Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 the figures. LC designed and co-supervised the project, and wrote the manuscript. CT co-supervised the project and assisted with array analysis and interpretation. AF provided assistance for running SES3D and data interpretation. LC was supported by a Vidi grant from The Netherlands Organisation for Scientific Research (NWO) on grant number 016.Vidi.171.022, as well as funding from the European Research Council under the European Union's Seventh Framework Programme 990 (FP/2007(FP/ -2013/ERC Grant Agreement no. 320639. FD was supported by the Ministry of Science and Technology (MOST) of Taiwan grant 107-2116-M-001-010. We thank Jeroen Ritsema and an anonymous reviewer for their comments, which helped to improve the manuscript significantly. Wookey, J., Stackhouse, S., Kendall, J., Brodholt, J. & Price, G.D., 2005. Efficacy of the post-perovskite phase as an explanation for lowermostmantle seismic properties, Nature, 438, 1004. Wysession, M.E., Lay, T., Revenaugh, J., Williams, Q., Garnero, E.J., Jeanloz, R. & Kellogg, L.H., 1998

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

A P P E N D I X I : D E TA I L S O F N U M E R I C A L S I M U L AT I O N S
The numerical set up used to model plumes in this study is essentially the same as that used in the global simulations of Deschamps et al. (2018). Conservation of mass, energy, momentum and composition are solved for a compressible, infinite Prandtl number fluid using the software StagYY (Tackley 2008). Details of the numerical techniques used to solve conservation equations may be found in Tackley (2008). The main properties of our models are detailed below, and parameter values are listed in Table A1. Fig. A1 shows the raw plumes, that is before adiabatic correction and re-scaling.

A1 Geometry and general physical properties
Conservation equations are solved in a spherical sector that covers an angular area a 45 • × 22.5 • and is sampled by 128 horizontal slices of 128 × 64 nodes each. The ratio between the radius of the core and the total radius is set to its Earth value, that isf = 0.55. The bottom and surface boundaries are free slip. The system is heated both from the bottom and from within, with a total mantle heating rate equivalent to a surface heat flux of 10 mW m -2 and a ratio of internal to basal heating close to 30 per cent. Compressibility generates additional sinks and sources of heat that are controlled by the dissipation number, Di, which varies with depth. We fixed the surface value of this number to Di S = 1.2, and the depth variations of thermal expansion imply that its volume average is equal to 0.43. Because the fluid properties (density, viscosity, thermal diffusivity, and thermal expansion) are allowed to vary throughout the system, the definition of the Rayleigh number is non-unique. In our simulations, we prescribed a reference Rayleigh number Ra 0 , defined at surface values of the thermodynamic parameters and reference viscosity η 0 . Here, Ra 0 is set to 3.0 × 10 8 , leading to an effective Rayleigh number (i.e. the Rayleigh number at the volume average viscosity) of about 10 6 for both TC1 and TC2.

A2 Thermochemical field
Thermochemical models distinguish two types of material, regular mantle and a chemically distinct, or primordial, material. The latter accounts for chemical heterogeneities that may be present at the bottom of the mantle as a result of early differentiation. The compositional field is modelled with a collection of 30 million of tracers, leading to an average number of tracers per cell of about 30, which is enough to properly model entrainment (Tackley & King 2003). Tracers are of two types, modelling the regular mantle and primordial material, respectively, and are advected following a 4th order Runge-Kutta method. At each time step, the compositional field is inferred from the concentration C of particles of primordial material in each cell, and varies between 0 for a cell filled with regular material only, and 1 for a cell filled with primordial material only. The exact nature of the compositional field is not prescribed a priori (except for its density excess) and can be fixed by the user during post-processing. Here, we assumed that the regular mantle is pyrolitic, and that the chemically distinct (primordial) material is enriched in iron oxide by 4 per cent (see main text). The primordial material is initially distributed in a basal layer. The thickness of this layer is controlled by the volume fraction of dense material, X prim , which we fixed to 3.5 per cent for both TC1 and TC2. The primordial material is assumed to be denser than the regular (pyrolitic) mantle, and the density contrast between the two materials is controlled by the buoyancy ratio, here defined with respect to a reference density that increases with depth following a thermodynamical model of Earth's mantle, where ρ c (z) is the density contrast between dense and regular material, α S the surface thermal expansion, ρ(z) the reference density at depth z, and T S the super-adiabatic temperature jump. The buoyancy ratio is set to B z = 0.23 in TC1, and B z = 0.18 in TC1. At the bottom of the system, and taking α S = 5.0 × 10 −5 K −1 , ρ bot = 4950 kg m −3 and T S = 2500 K, this leads to a density contrast between dense and regular material around 142 kg m −3 for TC1 and 110 kg m −3 for TC2.

A3 Viscosity
Viscosity is allowed to vary with depth, temperature, and composition. An additional viscosity ratio η 660 = 30 is added at the 660-km-phase transition. Furthermore, to avoid the formation of stagnant lid at the top of the system, we impose a yield stress. The viscosity is then fully described by where is the yield viscosity, and The yield viscosity (eq. A3) is defined from the yield stress, σ Y = σ 0 + σ z P, and the second invariant of the stress tensor,é. The Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 yield stress is set to σ 0 at the surface and increases with pressure following a gradient σ z . In eq. (A4) η 0 is a reference viscosity, H the Heaviside step function, z the depth, D the mantle thickness, T S the super-adiabatic temperature difference across the system, and T off the temperature offset, which is added to the temperature to reduce the viscosity jump across the top thermal boundary layer. The reference viscosity η 0 is defined for the surface value of the reference adiabat (i.e. T as = 0.64 T S ), and at regular composition (C = 0). The viscosity variations with temperature are controlled by E a , modelling the activation energy. To quantify the thermally induced increase of viscosity, we define a potential thermal viscosity ratio as η T = exp(E a ). However, due to the adiabatic increase of temperature and to the temperature offset, which we fixed to T off = 0.88 T S , the effective top-to-bottom thermal viscosity ratio is smaller than η T by about two orders of magnitude. Here, we fixed E a to 20.723, corresponding to η T = 10 9 . The viscosity variations with depth are controlled by V a , modelling the activation volume. Here, we fixed the value of V a to 2.303, leading to an increase of viscosity from top to bottom by a factor 300. Note that this increase includes the viscosity jump at 660 km, but excludes the decrease due to adiabatic increase of temperature and the thermally-induced increase in thermal boundary layers. The viscosity variations with composition are controlled by the parameter K a , and the viscosity ratio between primordial and regular material (or chemical viscosity ratio) is given by η C = exp(K a ). In this study, we impose primordial material to be more viscous than regular material with η C = 30, accounting for the fact that if dense material is enriched in bridgmanite (Trampert et al. 2004;Mosca et al. 2012), it may be more viscous than surrounding mantle (Yamazaki & Karato 2001).

A4 Phase changes
The transformation of ringwoodite into bridgmanite and ferropericlase at 660 km is modelled with a discontinuous phase transition controlled by defining a point on the phase boundary and a Clapeyron slope, 660 . Here, we imposed z = 660 km and T = 1900 K as anchor point, and 660 = −2.5 MPa K -1 . The density contrast at the phase transition is ρ 660 = 400 kg m -3 and is scaled with the surface density. Lateral deviations in the transition depth are determined using the phase function approach of Christensen & Yuen (1985). Combined with the 660-km viscosity increase (from upper to lower mantle), the 660-km phase change has a strong influence on the geometry of the plumes. This transition acts as a negatively buoyant barrier, which results in a spreading of the plume conduit beneath this boundary, and a thinning above it. The phase transition to post-perovskite, which is expected to happen in cold regions in the deep mantle, is not included here. Because this phase change has only a limited effect on the generation of plumes (e.g. Figure A0. Snapshot of non-dimensional, uncompressed temperature field from numerical simulations for (a) TC1 and (b) TC2. Isosurface value is = 0.55 in both case. Note that calculations of convection account for the adiabatic correction, but that output thermal distributions do not. For application to plume detection, we added again this contribution to the output temperature field, and rescaled this field with the super-adiabatic temperature jump following eq. (A5). Li et al. 2015), this omission may not alter our plume simulations and the seismic detection of these plumes.

A5 From output to real temperatures
All calculations are done in non-dimensional units. Adiabatic effects on temperature are taken into account when solving the energy and momentum conservation equations, but for practical reasons, output temperature fields do not include these effects. For application to plume detection in Earth mantle, output temperature must be rescaled and corrected with the adiabatic increase of temperature with pressure. The 'real' temperature at a given location, T(x,y,z), is then obtained from the non-dimensional, uncompressed temperature,T (x, y, z), following: where T top is the surface non-dimensional temperature, here fixed to 0.12, which is equivalent to a dimensional surface temperature of T surf = 300 K, T S = 2500 K the superadiabatic temperature jump, and a(z) the adiabatic correction at depth z. This correction is given by where Di S is, again, the surface dissipation number, and α(z) and C P (z) the thermal expansion and heat capacity as a function of depth. These two functions are defined as part of a reference thermodynamic model involved in the compressible form of conservation equations (Tackley 1998). Practically, α decreases by a factor 5 from the surface to the core-mantle boundary (CMB), while C P is constant with depth. The adiabatic correction defined in eq. (A6) then varies from 1.0 at the surface to about 1.55 at the CMB.

APPENDIX II: SUPPLEMENTARY FIGURES
This section contains the Figures A0 to A7.
Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Figure A1. Example horizontal temperature slices through TC1 and TC2. These slices are used to estimate the width of the plume and its thermal anomaly, as recorded in Table 1. Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Figure A2. Example of a 'bad' reference model: Percentage difference of Vp, Vs and density for TC1 relative to the reference, where in this case the reference was a 1-D vertical profile at the edge of the model (i.e. far from the plume centre). The seismic signal is dominated by differences in the transition zone discontinuities rather than the plume structure.
Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020 Figure A3. wave speed and density perturbations of the plume model TC2 from two different 1D reference models: AK135 (top row) and adiabatic pyrolite (bottom row).The left column shows the perturbations along a 1-D vertical profile far from the plume centre. The right column shows the perturbations along a 1-D vertical profile through the centre of the plume. With an adiabatic pyrolite reference, the profiles are dominated by the plume structure rather than systematic offsets in the transition zone discontinuities. Note that the apparent double-discontinuity at the base of the model is an artefact, with the chemical pile in TC2 causing the upper discontinuity and a phase change from bridgmanite to post-perovskite in the reference model causing the lower discontinuity.
Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020   TC2 shows a strong arrival at ∼850 s with a low slowness, that is not present in the reference model. This arrival appears to come from reflections off the chemical pile at the base of the model. Downloaded from https://academic.oup.com/gji/article-abstract/219/Supplement_1/S167/5544366 by Universiteitsbibliotheek Utrecht user on 03 February 2020