Temporal variations in plume flux: characterizing pulsations from tilted plume conduits in a rheologically complex mantle

SUMMARY

1 INTRODUCTION

Contrary to expectations of steady upwelling in the stems of mantle plumes (e.g. Morgan 1971), estimates of magmatic output through time along age-progressive hotspot volcano chains indicate variations of more than an order of magnitude (e.g. Wessel 2016). These variations can be smooth such as along the Louisville hotspot chain which, after a period of relatively constant magmatic output from ∼80 to ∼40 Ma, exhibited a steady decrease in extrusive volcanism (Beier et al. 2012; Vanderkluysen et al. 2014). Alternatively, variations can be punctuated like the rapid but temporary increases in magmatic production observed along the oldest portions of the Ninety East volcanic ridge (Coffin et al. 2002). Finally, estimates from more than a dozen hotspots including Hawaii, Iceland, Amsterdam-St Paul and Tristan da Cunha (O'Connor et al. 2000; Van Ark & Lin 2004; Vidal & Bonneville 2004; Maia et al. 2011; Wessel 2016; Morrow & Mittelstaedt 2021) suggest that variations commonly take the form of regular, likely periodic, pulses of magmatic output over millions of years. The periodicity of these pulses provides clues to the yet unconstrained mechanisms driving variable melt production within plume stems.

Studies that assess the periodicity of hotspot magmatic output often resolve strong spectral signals, commonly at periods on the order of a few to tens of millions of years. For example, along the Hawaiian hotspot chain estimates find volcanic production varies from 0.05 to 0.25 km³ yr⁻¹ at periods between ∼2 and 30.8 Myr (Van Ark & Lin 2004; Vidal & Bonneville 2004; Wessel 2016; Morrow & Mittelstaedt 2021). In the South Atlantic, volcanic production at the Tristan da Cunha and St Helena hotspots also shows periodic signals (Adam et al. 2005). Of these, Tristan da Cunha shows two primary periodicities, one between 10 and 20 Myr and another at ∼5 Myr. In contrast, the St Helena hotspot track shows two punctuated increases in volcanic production imprinted on top of a ∼5 Myr periodic variation. The observed range in periodicities of hotspot magmatic output suggests that multiple mechanisms likely influence temporal variability of plume conduits. Yet, the common ∼5 Myr period shared by several hotspots also points to a common mechanism for magmatic variability.

Mantle plumes that rise from the core–mantle boundary (CMB, Table 1) with a constant supply of material from the thermal boundary layer (TBL) and no external disturbances will form a steady, axisymmetric conduit (e.g. Albers & Christensen 1996). However, such ideal conditions are unlikely in the Earth's mantle and several processes can perturb the plume conduit, including solitary waves generated in the TBL (Schubert et al. 1989; Olson 1990), entrainment of compositionally different material by the plume (e.g. Iceland; Hanan & Schilling 1997) and plate-motion induced tilting of the plume conduit in the upper mantle (Whitehead 1982). Of particular interest here is the development of instabilities in tilted plume conduits. Laboratory experiments find that increasing conduit tilt in compositional plumes decreases plume vertical velocity and increases plume radius (Whitehead 1982). An increase in the radius of a sufficiently tilted plume (>∼60° from vertical) enhances the component of buoyancy oblique to the conduit, promoting instability formation and a consequent vertical pulse of material (Whitehead 1982). This seemingly large tilt may be applicable to plumes in the mantle. From a study of 44 possible hotspots, 3 separate models for calculating plume tilt find that multiple hotspots (e.g. Galapagos and St Helena) may have a tilt above 60° (Steinberger 2000). However, even above this critical tilt, scaling analyses based upon numerical and laboratory studies suggest that plume conduits may remain stable under mantle conditions (Kerr et al. 2008; Mériaux et al. 2011). Yet, previous studies examining instability growth in tilting plumes used one or multiple simplifying assumptions such as constant or only temperature-dependent rheology, a simplified plume conduit geometry and a uniform mantle structure (e.g. no phase changes), making extrapolation to mantle conditions difficult.

In addition to the potential role of complex rheology, the mantle transition zone can alter plume dynamics. As a plume rises through the mantle it interacts with multiple phase transitions, where mineral structures change in response to differences in pressure and temperature. Of these, the two most significant are the 660-km transition (Bridgmanite and Ferropericlase to Ringwoodite), an endothermic reaction with a Clapeyron slope of –4 to –0.4 MPa K^–1 (Ito et al. 1990; Ghosh et al. 2013; Faccenda & Dal Zilio 2017), and the 410-km transition (Wadsleyite to Olivine), an exothermic reaction with a Clapeyron slope between 1 and 4 MPa K^–1 (Akaogi et al. 1989; Katsura et al. 2004). Together these two transitions bound a region known as the transition zone (Flanagan & Shearer 1998), and work to inhibit (660-km) or promote (410-km) plume upwelling. The degree of inhibition or promotion depends on the magnitude of the Clapeyron slope (Schubert et al. 1975). In addition, a large jump (∼30×) in viscosity at approximately the depth of the 660 km phase change may anchor the lower portion of an upwelling plume and facilitate tilting of the conduit in the upper mantle.

Numerous studies have examined the dynamics of thermal plumes, as well as perturbations to those dynamics, including the role of the transition zone (Nakakuki et al. 1994; Schubert et al. 1995; Bossmann & Van Keken 2013), deep piles of anomalous material (Steinberger & Torsvik 2012) and the impact of grain size evolution (Dannberg et al. 2017). However, no consensus exists as to the processes responsible for variability in plume upwelling. In this study, we aim to isolate the role of the mantle transition zone and plate motion on instability formation in tilting plumes as one potential mechanism for explaining observed variations in volcanic production along hotspot tracks (e.g. Kerguelen, Hawaii and Tristan da Cunha). In contrast to previous studies, we simulate plume evolution in a mantle with a complex rheology (diffusion and dislocation creep) and phase transitions. Our model results indicate that when plate motion-driven mantle shear tilts plume conduits to >∼67°, the conduit destabilizes, initiating periodic pulsing (∼2–6 Myr periods) and consequent variability in plume melting. The periodicity is not a strong function of the Clapeyron slopes of the 410-km phase transition but depends strongly on the viscosity ratio between the ambient mantle and plume conduit, the overlying plate velocity, and the effectiveness of dislocation creep.

2 METHODS

2.1 Governing equations

2.2 Phase function

2.3 Reference profiles

2.4 Rheology

We include temperature jumps associated with latent heat release during phase changes in the temperature reference profile. The temperature jumps depend on the Clapeyron slope, resulting in small differences in reference temperature profiles between models. To account for this in our viscosity profile, the diffusion creep constant pre-factor was varied to account for changes in temperature related to different 410-km Clapeyron slopes (Table S1), but activation energy and activation volume were unchanged. Pre-factor values were chosen to bring viscosity profiles as close as possible to a reference profile, especially in the upper mantle (the focus of this work), but profiles still differ slightly from each other. Fig. 1 shows viscosity, temperature, density, thermal expansivity and thermal conductivity profiles at time step 0 covering the range seen within all the models due to the temperature jumps.

Figure 1.

Model parameter profiles at time step 0, including (a) viscosity (linear rheology), (b) temperature, (c) density, (d) thermal conductivity and (e) thermal expansivity. The solid line shows the coldest run (410-km = 4 MPa K^–1) and the dashed line the warmest run (410-km = 1 MPa K^–1). Note that thermally driven density changes are governed by deviations in temperature from the adiabatic temperature profile which includes temperature jumps at phase changes. Thus, the density profiles at time step 0 are similar regardless of differences in initial temperature profiles.

Open in new tabDownload slide

2.5 Model setup

To numerically simulate an isolated thermal plume, we use a 2-D shell geometry with an opening angle of θ = 90° and mantle thickness of 2855 km (Fig. 2). As the focus of this work is on the effects of phase transitions and plate shear on an upwelling plume, we do not include an upper thermal boundary layer and set the upper model surface to a depth of 35 km. The choice of 35 km is equivalent to a ∼10 Ma oceanic lithosphere across the model surface, but the exact choice of this depth is unlikely to significantly change the dynamics of our model (see Section 4.5). Therefore, along our upper surface, we impose a surface pressure of 1.18 × 10⁹ Pa, equivalent to a 35-km-thick lithosphere with a density of 3300 kg m⁻³ and a constant temperature boundary condition with T(surface) = 1600 K, equivalent to expected mantle potential temperatures (Courtier et al. 2007; Herzberg et al. 2007). Thus, we ignore the role of a constantly increasing plate thickness and avoid possible complications due to foundering of destabilized lithosphere.

Figure 2.

Initial model setup. Colours indicate initial model temperature (blue to red). Physical boundary conditions include a free-slip lower boundary, fixed lateral velocity upper boundary and mixed conditions on the vertical boundaries. At depths below 485 km, the vertical boundaries are free slip while for shallower depths inflow is specified on boundary θ = 90°, and outflow on the boundary θ = 0° (arrows and inset in upper right). Along the outer boundary we impose a constant boundary-parallel velocity of between 2 and 10 cm yr⁻¹ to simulate plate motion. White lines indicate the phase transitions at depths of 410 and 660 km.

Open in new tabDownload slide

Mechanically, we impose a constant, uniform, surface-parallel velocity between 2 and 10 cm yr⁻¹, a free-slip inner (bottom) boundary and mixed boundary conditions on the left and right edges of the model. On the left and right boundaries, a free-slip condition is set below 450 km depth (see Fig. S1 for differences based on this chosen depth) and a linearly decreasing normal velocity (i.e. Couette flow) is used to prescribe inward (left boundary) and outward (right boundary) flow (e.g. Turcotte & Schubert 2013) at shallower depths. The inflow and outflow fluxes balance and account for the flow of upper mantle material far from the plume due to the lateral surface velocity. Thermally, the surface and base of the mantle are isothermal (T_top = 1600 K, T_bot = 3200 K) and the sides have zero lateral temperature gradient (no heat flux).

Initial conditions for the model include an adiabatic temperature gradient, Clapeyron slope dependent jumps in temperature (up to 81 K at the 410 km transition and –42 K at the 660 km transition) and density (165 and 330 kg m⁻³) due to the 410 and 660 km phase transitions, a lower thermal boundary layer (ΔT = 570 K) that has evolved for 300 Myr, and a localized 600 km wide thermal perturbation centred on the bottom boundary (ΔT = 550 K) to initiate plume upwelling (Fig. 2). The 410 and 660 km phase transitions are included at model depths of 375 and 625 km to account for the imposed 35 km depth of the model upper surface. Hereafter all depths in the text will assume that the upper surface depth is 35 km (e.g. the depth to the phase transitions will be referred to a 410 and 660 km). Phase transitions thicknesses are set to 30 km, and Clapeyron slopes are varied between 1 and 4 MPa K^–1 for the 410 km and set at a fixed value of –1 MPa K^–1 for the 660 km transition.

Models were run on 64 processors for ∼100 hr (until reaching a stopping criteria of 400 Myr model time). Model resolution varies spatially initially and adaptively throughout each simulation. The coarsest elements measure ∼90 km across and the more refined elements ∼6 km across. Adaptive refinement of the mesh was based on viscosity and non-adiabatic temperature (ΔT), yielding the finest resolution within and around the plume conduit. For phase transitions, we imposed a maximum element size of ∼24 km to avoid artefacts in latent heat produced by flow. As we are interested in plume dynamics in the upper mantle, the most refined region of the model is between 32° and 58°, and >700-km model depths. Solver tolerances were set as 10⁻⁸, 10⁻⁵ and 10⁻¹² for the linear, non-linear and temperature solvers, respectively.

2.6 Model calculations

2.6.1 Plume volume flux

To calculate the volume flux of the plume in the upper mantle we assume that any outward radial flow near the plume is primarily due to plume upwelling. First, we interpolate model velocities and non-adiabatic temperatures to a depth contour of 175 km below the model surface (i.e. a mantle depth of 210 km including our assumed lithosphere thickness). Secondly, we find the maximum non-adiabatic temperature along this depth contour and assume this is the plume location. Finally, we assume that any outward flow within ±1° in θ of this maximum temperature and along this depth contour is related to upwelling plume and integrate the radially outward (i.e. locally vertical) velocities within range (∼215 km arc length at this depth).

2.6.2 Plume periodicity

We assess the periodicity of upper mantle pulsing within the model by examining the upper mantle plume flux and utilizing a Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) to identify spectral peaks in our plume flux time series with a statistical significance >99 per cent. We apply this method in two separate ways, (1) the entire volume flux data from Upper Mantle Pulsing (UMP) initiation to the end of the model is used to get an ‘overall’ periodicity. (2) To examine the temporal evolution of we use the same time frame, but instead calculate periodicity within a 25 Myr wide Running Window (RW) that is centred at 10 Myr intervals. In this method the first interval is chosen so that no volume flux data before UMP initiation is included. The final interval is chosen so that the 25 Myr window does not exceed the model end time.

2.6.3 Plume tilt

To quantify changes in plume stem tilt from vertical we examine the maximum plume tilt averaged over 25 km increments across a depth interval of 75–575 km below the model surface (i.e. mantle depths of 110–610 km including our assumed lithosphere thickness). By limiting our depth to a minimum of 70 km, we avoid problems of extremely high tilts where the plume is advected horizontally along the top model boundary. To find the plume tilt we first examine depth contours at 5 km intervals and extract the location of the maximum temperature in radial coordinates under the assumption that it represents the plume. We note that because the mesh refines adaptively, not all depth contours will have a data point and the radial distance between points may differ. Additionally, to avoid complications that may arise from highly tilted conduits due to this method (e.g. maximum temperatures along a depth contour that are not associated with the plume), we limit the maximum tilt to ∼84°. This does not appear to affect the results presented here as the tilts generally fall below ∼75°. Using these maximum temperatures, we skeletonize the plume down to a line of data points (e.g. Fig. S2). From here, we fit a 10th degree polynomial to the line with data points at 5 km intervals. Next, we determine the tilt value between each 5 km depth increment. We then run a five-point averaging filter along the slope values to calculate an average tilt value at each depth increment. Finally, we use the maximum value to describe the tilt of the plume stem.

2.6.4 Plume temperature, viscosity and mantle viscosity

To evaluate how the plume evolves through time, we track the plume temperature, plume viscosity and mantle viscosity. First, we interpolate the model data to find the temperature and viscosity along the 340-km depth contour beneath the model surface (i.e. a mantle depth of 375 km including our assumed lithosphere thickness). This depth contour is chosen as it is generally deeper than where plume pulsing occurs ensuring that there is a constant supply of plume material to evaluate. Next, we find the maximum plume temperature along this depth contour and similarly use this point to describe the plume viscosity. For the mantle viscosity and temperature, along the same depth contour we average the viscosity from ∼200 to 1200 km upstream in theta from the plume maximum temperature.

2.6.5 Melt production of a mantle plume pulse

To evaluate the likely surface expression of plume pulses, we integrate the total deviation in modelled melt production caused by each mantle plume pulse. To simulate melting, we assume that the fraction of material undergoing melting in our models is pressure and temperature dependent and follows the empirically determined anhydrous melting equations of Katz et al. (2003). Numerically, ASPECT tracks the melt fraction F as a compositional field that is output as a globally integrated statistic every time step. To determine the volume of melt produced in each time step we first smooth this data using a running average over 10 time steps (∼100 to 150 kyr). Next, we calculate the difference in F from the previous time step and spatially integrate this difference across the melt compositional field. To determine the melt production of a single plume pulse, we isolate the temporal limits of a single plume pulse by locating local minimums in melt production (Fig. S3). Using a straight line fit between each pair of minimums, we define melt production by a pulse as the melt production values above that line. The area above the line is integrated, yielding the total melt produced by the pulse in km³ km⁻¹. Average pulse melt production rate is this integrated value divided by the time elapsed between the two minimums.

3 RESULTS

We present results from 18 model simulations where we vary parameters that may govern periodic variations in upwelling of tilted plume conduits including mantle rheology, the Clapeyron Slopes of the 410 km depth mantle phase transitions, overlying plate speed, and plume excess temperature. We compare all simulations to a reference model that consists of a composite rheology (eq.. 18), 410-km phase transition Clapeyron slope of 3 MPa K⁻¹, overlying plate speed of 8 cm yr⁻¹ and an initial plume temperature anomaly of 550 K and a CMB temperature anomaly of 570 K. First, we examine the role of mantle rheology through a set of six simulations. These simulations consist of four model runs that have a linear rheology (i.e. eq. 16) where we vary the Clapeyron slope of the 410-km transition from 1 to 4 MPa K⁻¹, and two model runs that use a composite rheology and reduce or increase the dislocation pre-factor of the upper mantle by a factor of 10. Next, we examine the role of overlying plate speed by comparing additional simulations using values of 2, 4, 6 and 10 cm yr⁻¹ to the reference model. Finally, we examine the role of initial plume temperature by modifying the initial plume and CMB temperature anomaly of the reference model to 450, 650 or 750 K. Model outputs are evaluated on their general evolution, vertical plume flux in the upper mantle, types of plume pulse initiation, periodicities of plume pulses and changes in melt production associated with pulses.

Note that in our results and discussion below we refer to two types of model time. First, we cite standard times, where the zero corresponds to the start of a simulation; we note these times as model time (MT). Additionally, we are interested in the timescale for development of tilting and instability formation in our models as they relate to observed ages of hotspot islands so we also cite times after plume arrival (APA) in the upper mantle.

3.1 Reference model

To facilitate description of varying model cases, we first describe the evolution of a reference case with composite rheology, a plate speed of 8 cm yr⁻¹, an initial plume anomaly of ΔT 550 K, and Clapeyron slopes of the 410 and 660 km phase transitions set to 3 and –1 MPa K⁻¹, respectively (Fig. 3). We then present results from simulations with different Clapeyron slopes, overlying plate velocities, initial plume temperature anomalies and upper mantle dislocation pre-factors by comparison to this reference model.

Figure 3.

Simulation results for the reference model. The simulation shows variability in vertical plume volume flux (a) associated with destabilization of plume conduits (b) that tilt >67° from the vertical. Statistically significant periods (99.9 per cent statistical significance) of plume flux variability over a running window (considered volume flux for the red dot is indicated by the red arrows in a) indicate a general increase in periodicity through time (c). Melt produced by discrete volume flux pulses (d) are distinguished by alternating shades of blue. (e) Overall periodicity taken over the entire time of UMP (dashed line, α = 0.01 indicating 99.9 per cent statistical significance). (f–j) Model snapshots at different times (indicated by green triangles in a and b), where the contour for volume flux calculations is shown in purple. Contour for plume parameters is shown in blue. White contours show the non-adiabatic temperature from 25 to 200 K, at 25 K intervals. Arrows show the direction of the velocity field (not scaled by magnitude).

Open in new tabDownload slide

In the reference case (Video S1), the plume head forms and begins to rise from the core–mantle boundary (CMB) at ∼90 Myr MT. The plume head increases in speed as it rises through the lower mantle, reaching a maximum velocity of 11.5 cm yr⁻¹ before reaching the 660 km phase transition at ∼139 Myr MT. Upon reaching the base of the transition zone, the plume head slows as it crosses the endothermic 660 km depth phase change. Continuing across the exothermic 410 km depth phase change, the plume radius decreases (∼190 to ∼42 km when examining the non-adiabatic temperature (>50 K) at 142 Myr MT and depths of 625 and 375 km, respectively), and the maximum velocity briefly increases to 207.8 cm yr⁻¹. The plume head rapidly crosses the upper mantle (∼2 Myr) and, initially, spreads symmetrically beneath the model surface. Shear due to the imposed lateral velocity along the upper model surface (i.e. plate motion) causes plume material to subsequently spread in the direction of plate motion. Over the next ∼28 Myr (142 Myr to 170 MT), as the remainder of the plume head crosses the transition zone, the plume conduit responds to plate shear by tilting from vertical to a maximum angle of 72° (Fig. 3b). Concurrently, plume temperature and volume flux increase in the portion of plume stem beneath the 660 km transition. The increase in volume flux travels through the tilted plume stem and results in three large pulses rising from the 660 km transition. Each large pulse reorients the plume stem closer to vertical in the uppermost mantle.

Hereafter, we refer to the generation of these deep pulses as Transition-Zone Pulsing (TZP, Video S1, ∼0:30).

Following TZP (∼173 to 193 Myr MT), the plume stem tilts again, reaching a maximum tilt of 69° by 201 Myr. After this, an instability develops within the uppermost mantle that leads to a pulsing of the plume stem. (∼202 Myr MT, Video S1, ∼0:39). To highlight the difference to TZP, we refer to this shallower destabilization process as upper mantle pulsing (UMP). Destabilization of the plume conduit as described here for the reference case is the first of two mechanisms we observe for initiating UMP in our models and is referred to as UMP type 1 (Figs 4a–d). UMP continues for the remainder of the model run. Initially, pulses form at a model depth of ∼150 km, but deepen over time, reaching a maximum model depth of ∼310 km by the end of the reference case.

Figure 4.

Example of UMP initiation for both types. (a–d) UMP type 1 initiation by shearing of the plume stem until a high enough tilt is reached for instabilities to form (model run with composite rheology, and a 410-km Clapeyron slope of 1 MPa K^–1; Video S4). (e–h) UMP type 2 initiation by destabilization from deeper TZP pulses (model run with composite rheology, a 410-km Clapeyron slope of 3 MPa K^–1, and an overlying plate velocity of 6 cm yr^–1; Video S6).

Open in new tabDownload slide

Next, we examine periodicity and melt production related to UMP in the reference model. Periodicity is examined in two ways as described in Section 2.6.2 that we refer to as running window (RW) and overall periodicity. We find that when examining the overall periodicity from the start of UMP to model end, two signals appear at 4.0 and 3.7 Myr (Fig. 3e, Table S2), with the most significant signal at 4.0 Myr. By instead looking at RW periodicities, we see a gradual increase in periodicity with time ranging from a minimum of 2.7 Myr at the start of UMP to a maximum of 4.1 Myr, with periodicity remaining relatively constant at ∼4.0 Myr from ∼150 to 225 Myr APA. The lowest periodicities occur near UMP initiation, and the maximum periodicities at the end of the model run (Fig. 3c). For the total melt produced by a UMP pulse, we find ranges between 217.2 and 845.1 km³ km⁻¹, with a median of 508.3 km³ km⁻¹ (Table S3).

3.2 Effects of rheology on plume evolution

In the reference model we use a composite rheology with an uppermost mantle dislocation pre-factor of 8.33 × 10⁻¹⁷ Pa⁻ⁿ s⁻¹. Here we discuss the effects of model rheology by comparing six simulations to the reference model. Four of these simulations utilize solely diffusion creep with no dislocation component wherein we vary the 410-km Clapeyron slope from 1 MPa K⁻¹ to 4 MPa K⁻¹, one simulation that has an uppermost mantle pre-factor reduced by a factor of 10 (8.33 × 10⁻¹⁸), and one simulation where we increase the dislocation pre-factor by a factor of 10 (8.33 × 10⁻¹⁶).

The four simulations that utilize solely diffusion creep evolve similarly to the reference run, exhibiting TZP and UMP, but differ in the time it takes the plume to reach the upper mantle. In all cases, the plume takes longer to reach the upper mantle with increasing Clapeyron slope (Table 4). The initiation of UMP also differs between the runs. While the simulation with a Clapeyron slope of 1 MPa K⁻¹ initiates by UMP type 1 like the reference model, in the other cases UMP initiation is connected to instabilities during TZP. Smaller TZP pulses rise and greatly increase plume tilt in the uppermost mantle. This leads to large UMP pulses that eventually stabilize into continuous UMP. We refer to UMP initiation by TZP as UMP type 2 (Figs 4e–h; Video S2, ∼0:35). The two simulations with a 10 times reduction or increase to the dislocation pre-factor (Video S3) show similar timings to the reference model, initiate by UMP type 1 (Table 4), and exhibit TZP and UMP.

Next, we investigate the role of non-linear rheology on UMP periodicity within the simulations. All simulations show a general increase in RW periodicity over time, with values across all models ranging between 1.6 and 6.7 Myr (Table S2), similar to the reference model (Fig. 5a). When examining the overall periodicity signals, the models show multiple signals clustered between 3.1 and 5.2 Myr. Comparing the linear rheology model most similar to the reference model (i.e. with a 410-km Clapeyron slope of 3 MPa K⁻¹) to the other model runs with varying dislocation pre-factor, we see fewer signals as dislocation becomes more effective (e.g. 6 signals with linear rheology and 1 signal with a pre-factor 10 times larger than the reference model; Table S2). Additionally, if we compare only the most significant signal from each model, there is a trend of decreasing periodicity as dislocation creep becomes more effective (Fig. 6a).

Figure 5.

Periodicity versus time (APA) for Category 2 simulations with different (a) dislocation pre-factor, (b) 410-km Clapeyron slope, (c) overlying plate velocity and (d) plume/CMB temperature anomalies. Lines connect the most significant periodicity for each time. In all plots, the black line represents the reference model.

Open in new tabDownload slide

Figure 6.

Pulsing plume periodicity against dislocation pre-factor (a), 410-km Clapeyron slope (b), overlying plate velocity (d) and Initial thermal perturbation (d). Average melt produced per pulse against dislocation pre-factor (e), 410-km Clapeyron slope (f), overlying plate velocity (g) and Initial thermal perturbation (h). Black lines indicate composite rheology, light red lines indicate linear rheology.

Open in new tabDownload slide

Finally, we examine the effect of rheology on melt production. As most of our models consider a composite rheology, and dislocation creep is thought to be important in the uppermost mantle (e.g. Karato & Wu 1993; van Hunen et al. 2005), we only calculate the melt for composite rheology models. For the models with increasing dislocation effectiveness, pulse sizes vary between a total range of 46.6–845.1 km³ km⁻¹ (Table S3). The median of pulses for each run suggests no clear trend with dislocation creep effectiveness (Fig. 6e), with pulse sizes of 453.9, 508.3 and 184.0 km³ km⁻¹ for the reduced, reference and increased dislocation pre-factor runs, respectively.

3.3 Effects of 410-km Clapeyron slope on plume evolution

We ran three simulations varying the 410-km Clapeyron slope from 1 to 4 MPa K⁻¹ to test how plume evolution, UMP periodicity, and UMP melt production are affected by the 410-km phase transition. These models evolve similarly to the reference model exhibiting both TZP and UMP, although timings of plume arrival and UMP initiation type differ between them. Like the linear rheology models, we see a trend of increasing plume arrival time with greater 410-km Clapeyron slope (Table 4). The UMP initiation type shows no trend but varies between the runs, with Clapeyron slopes of 1 MPa K⁻¹ (Video S4) and 3 MPa K⁻¹ (reference model) initiating by UMP type 1. Conversely, models with 410-km Clapeyron slopes of 2 and 4 MPa K⁻¹ initiate by UMP type 2.

Like the reference model, RW periodicity across all models varying Clapeyron slope general increase with time ranging from 2.4 to 4.3 Myr (Fig. 5b; Table S2). The overall periodicity shows multiple signals for each model that range from 2.0 to 4.1 Myr. When examining the most significant signal in the overall periodicity, there is a small increase with increasing Clapeyron slope, with a signal of 4.0 Myr for models with a Clapeyron slope of 1, 2 and 3 MPa K⁻¹, and signal of 4.1 Myr in the model with a Clapeyron slope of 4 MPa K⁻¹ (Fig. 6b).

Like pulse periodicities, we also found that pulse melt production was affected by the Clapeyron slope. The per-pulse melt production for all four simulations varies between 80.3 and 880.8 km³ km⁻¹ (Table S3). The median pulse for each run was 331.4, 437.0, 508.3 and 579.7 km³ km⁻¹ for Clapeyron slopes of 1, 2, 3 and 4 MPa K⁻¹, respectively, suggesting a general increase in pulse size with greater Clapeyron slopes (Fig. 6f).

3.4 Effects of overlying plate velocity on plume evolution

Although the reference model has an overlying plate speed of 8 cm yr⁻¹, plate speed is likely to be a key parameter controlling plume tilt which previous studies demonstrate is important for generating instabilities within the plume stem (e.g. Whitehead 1982). We ran five additional models varying the plate speed from 0 to 10 cm yr⁻¹ to quantify its effect on model evolution. In all models, plume arrival occurs at ∼140 Myr. In the case with no plate motion, no TZP or UMP occurs and the plume remains stable throughout the run. At 2 cm yr⁻¹, we observe a single, small TZP pulse, but UMP does not initiate and the plume remains stable (Video S5). The models with velocities of 4, 6 and 10 cm yr⁻¹ behave similarly to the reference model and exhibit TZP and UMP. For the 4 and 10 cm yr⁻¹ models, UMP is initiated by type 1 (as in the reference model). In the 6 cm yr⁻¹ model UMP is initiated by type 2 (e.g. Figs 4e–h; Video S6; Table 4).

Examining UMP periodicities, we find that changing overlying plate speed alters the simulated pulsing periods. Like the previous models, RW periodicities increase with model run time ranging from 2.6 to 6.2 Myr (Fig. 5c; Table S2). Each model shows multiple signals in the overall periodicity between 3.6 and 6.1 Myr. When examining only the most significant of these signals for each model, we see an increase in periodicity with slower plate speeds, with periodicities of 3.6 Myr at 10 cm yr⁻¹, 4.0 Myr at 8 cm yr⁻¹, 4.8 Myr at 6 cm yr⁻¹ and 6.1 Myr at 4 cm yr⁻¹ (Fig. 6c). Thus, increasing plate speed reduces the period of UMP in our models.

The volume of melt produced by UMP pulses is also likely affected by the overlying plate speed. Overall, individual pulses for the models vary within a large range, producing between 90.6 and 2342.7 km³ km⁻¹ of melt (Table S3). The median pulse for each run greatly increases with plate speed (Fig. 6g), from a median of 142.5 km³ km⁻¹ at a plate speed of 4 cm yr⁻¹, 195.7 km³ km⁻¹ at 6 cm yr⁻¹, 508.3 km³ km⁻¹ at 8 cm yr⁻¹ and 877.7 km³ km⁻¹ at 10 cm yr⁻¹.

3.5 Effects of initial plume temperature anomaly on plume evolution

Lastly, to understand how plume temperature affects model evolution we ran three additional simulations where we varied the initial ΔT of the plume and CMB from 450 K, 550 K (reference model), 650 K and 750 K (Video S7). In all cases, plume dynamics were similar, and each model exhibited TZP and UMP. Differences occurred in plume arrival time, with greater temperatures resulting in earlier plume arrival (e.g. arrival at 182.0 Myr with the coldest plume versus 90.5 Myr with the warmest plume) and UMP initiation type. Models with initial plume temperature anomalies of 550 and 650 K initiated by UMP type 1, and models with anomalies of 450 and 750 K initiated by UMP type 2 (Table 4).

The initial plume temperature anomaly affected UMP periodicities. Similar to previous models, RW periodicities increase with time ranging from 1.8 to 4.5 Myr across all models (Fig. 5d). Examining the overall periodicity shows multiple signals for each model ranging between 3 and 4.3 Myr (Table S2). The most significant signal for each model decreases from a periodicity of 4.3 Myr with an initial temperature anomaly of 450 K to 4.0 Myr at 550 K, 3.7 Myr at 650 K and 3.4 Myr at 750 K, suggesting that lower temperatures result in slower UMP pulsing (Fig. 6d).

Like the periodicities, changes in initial plume temperature anomaly affected the melt produced by a plume pulse. All models show melt production from a pulse that ranges in size from 49.9 to 940.0 km³ km⁻¹ (Table S3). However, the median pulse for each model decreases with greater initial plume temperatures (Fig. 6h), with median pulses of 550.5, 508.3, 396.8 and 349.9 km³ km⁻¹ for temperature anomalies of 450, 550, 650 and 750 K. This decrease in the melt produced in a given pulse with increasing temperature likely relates to the periodicity. The higher temperature models have a lower periodicity, suggesting a higher frequency of smaller pulses.

4 DISCUSSION

4.1 Initiation of upper mantle pulsing

In our simulations, UMP initiates in two ways: Type (1) the plume gradually tilts until instabilities form at the surface or Type (2) deeper pulsing from the transition zone continues in the upper mantle, destabilizes the tilted plume stem and initiates shallow pulsing. In both cases, the plume stem tilt smoothed using a 5 Myr running average exceeds and remains above a critical value while UMP is active. To isolate the relationship between tilt and instability development from other possible perturbations, we restrict the following discussion to UMP Type 1.

Plume tilt in UMP type 1 cases can generally be divided into two stages: (1) a gradual increase of tilt after TZP and prior to UMP initiation as the plume is sheared by overlying plate motion and (2) oscillations of the tilt with an amplitude of ∼10–20° during active UMP (Fig. 3b). During the first stage, the maximum tilt just prior to UMP initiation reaches values of 68–72°, with 5 of 8 cases tilting to >70° (Table 5). The only models that did not initiate UMP are those with plate velocities of 0 and 2 cm yr⁻¹. In these simulations the maximum tilts reached were 38° at 0 cm yr⁻¹ and 63° at 2 cm yr⁻¹ (we note that in the 2 cm yr⁻¹ case, the highest tilt within the model was 67° and occurred before a TZP pulse). Based upon these results, we suggest that a plume tilt >67° is needed for UMP to initiate, similar to the 60° proposed by Whitehead (1982) based upon analogue experiments.

Although our models predict that the initiation of UMP requires a minimum plume tilt >67°, after initiation the average tilt of the plume stem for all models that exhibited UMP ranges between 60° and 67° (with an average of 64°). In comparison, when examining a similar time frame (e.g. 80 Myr APA to the end of the model run) in the two models without UMP, they exhibited average tilts of 24° and 59° for the 0 and 2 cm yr⁻¹ cases, respectively. That UMP can continue at tilts less than the 67° tilt needed to initiate UMP suggests that it is easier to sustain UMP than to initiate it, though from the models presented here the minimum tilt needed to sustain UMP is unclear.

In addition to reaching a critical plume tilt above 67°, Type 1 UMP instabilities require 1–21 Myr to form (i.e. after reaching the critical tilt; Table 5), with a mean of 9 Myr. To conclude, we find that in the absence of external perturbations UMP initiates in plume stems that remain at tilts >67° for at least 1 Myr.

4.2 Controls on the periodicity of upper mantle pulsing

To investigate the controls on UMP periodicity, we restrict our discussion to model cases with a composite rheology as this includes the majority of our simulations and dislocation creep is believed to be important in the upper mantle (e.g. Karato & Wu 1993; van Hunen et al. 2005). In these models, UMP period increases with increasing model run time (Fig. 5), decreasing plume excess temperature (Figs 5d and 7a), and increasing plate velocity (Fig. 5c). Changes in the above parameters alter the temperature or strain rate in the upper mantle, the plume, or both. Over the course of a model run the temperature contrast between the plume and the ambient upper mantle decreases as the plume cools and/or the upper mantle is warmed by spreading plume material. Similarly, increasing the initial plume anomalous temperature increases the temperature contrast between plume and upper mantle. Additionally, faster plate velocities also increase the plume-mantle temperature contrast by accelerating the removal of anomalously warm plume material from the upper mantle and out of the model domain, maintaining a cooler upper mantle. Finally, in addition to cooling the upper mantle, faster plate motions will yield higher strain rates across the shearing depth of the upper mantle. The above changes not only change the temperature contrast between plume and mantle but will consequently alter the viscosity contrast between plume and mantle. Together, the above relationships suggest that the viscosity ratio (VR) between mantle and plume exerts a strong control on instability growth in the plume stem (Fig. 7), consistent with the conclusions of previous studies (e.g. Kerr et al. 2008; Mériaux et al. 2011).

Figure 7.

Plume pulsing periodicity of composite rheology models with varying Clapeyron slope, plate speed, or initial temperature anomaly. (a) Plume excess temperature with colours indicating viscosity ratio, with no applied scaling. (b) Viscosity ratio with colours indicating plume excess temperature, with no applied scaling. (c) Viscosity ratio with colours indicating plume excess temperature, with scaling applied to account for viscosity differences related to plate speeds (see Section 4.2) and a viscosity ratio cut-off below 6. (d) Viscosity ratio with colours indicating plume excess temperature, with scaling applied, using an exponential fit.

Open in new tabDownload slide

As the viscosity ratio is the primary control on the pulsing period, it is important to understand how the separate modelled parameters affect the ratio. Fig. 8 suggests that in general, lower 410-km Clapeyron slope (Fig. 8b) or higher initial temperature anomalies (Fig. 8d) have the largest effect on the viscosity ratio, with both increasing the ratio. Higher plume anomaly temperatures will result in greater plume temperatures in the upper mantle and lower plume viscosities, resulting in a lower viscosity ratio. Clapeyron slope likely similarly relates to plume temperatures. Because of the slightly different viscosity profiles related to the 410 km transition temperature jumps, models with lower Clapeyron slopes reached the 660-km transition faster resulting in a warmer plume. Conversely, mantle rheology (Fig. 8a) and overlying plate velocity (Fig. 8c) do not have a large effect on the viscosity ratio.

Figure 8.

Viscosity ratio versus time (APA) for (a) dislocation pre-factor, (b) 410-km Clapeyron slope, (c) overlying plate velocity and (d) plume/CMB temperature anomalies. In all plots, the black line represents the reference model.

Open in new tabDownload slide

4.3 The impact of non-linear rheology on the distribution of upper mantle pulsing periods

To assess the importance of non-linear rheology, we ran 3 simulations similar to the reference model (410-km Clapeyron slope of 3 MPa K⁻¹, plate velocity of 8 cm yr⁻¹, and initial plume temperature anomaly of 550 K) where we either increased or decreased the dislocation creep pre-factor by 10× in the uppermost mantle, or assumed a solely linear rheology model (i.e. only diffusion creep). A larger pre-factor results in the harmonic average between dislocation and diffusion creep (eq. 18) favouring the strain-rate dependent rheology. In other words, given the same strain-rate and temperature the dislocation creep related viscosity will be lower and have a larger contribution to the effective viscosity (eq. 18). Similarly, a smaller pre-factor results in a smaller contribution to the effective viscosity.

The two primary results from varying the effectiveness of the non-linear rheology are: (1) periodicity variance increases with decreasing effectiveness of dislocation creep (Table S2) resulting in more signals considered significant and (2) increasing the effectiveness of dislocation creep reduces UMP periodicity (Figs 5a and 6a). Based upon these results, the inclusion of non-linear rheology appears to stabilize the pulsing period, and further increasing its effectiveness reduces the UMP period. Thus, if dislocation creep is less effective in the upper mantle than assumed in this study, we may expect longer more diverse pulsing periods. On the other hand, if dislocation creep is more effective, then we would expect shorter periods.

4.4 The role of the upper thermal boundary layer in plume dynamics

The upper thermal boundary layer on Earth will perturb the spread of underlying plume material and may alter pulse formation due to UMP (e.g. Ribe 1996; Ito et al. 2003). Effects on instability growth could be caused by differences in layer thickness at different plume locations, layer growth (i.e. the oceanic lithosphere), or layer stability. Larger boundary layer thicknesses will deepen the top of the melting column and decrease the total volume of melt produced in each pulse (but not the periods of melt pulses). Similarly, a growing thermal boundary layer (i.e. the oceanic lithosphere <∼80 Ma; Stein & Stein 1994) promotes plume (and melt) flow up slope and towards younger lithosphere (e.g. Mittelstaedt et al. 2011), but the slope of this boundary decreases with increasing seafloor age as does its effect on the dispersal of plume material (e.g. Ribe 1996; Ito et al. 2003). These theoretical studies demonstrate that the upslope flow of plume material depends on viscosity and plate motion rate. Finally, the thermal boundary layer in the oceans reaches a constant thickness by ∼80 Ma, which may be maintained by destabilizations of the thermal boundary layer as the lowermost portions of the lithosphere ‘drip’ into the mantle (e.g. Dumoulin et al. 2005). Such drips may trigger UMP in a tilted plume, similar to how TZP triggers UMP for Type 2 cases. Although inclusion of a thermal boundary layer may impact the development of UMP, the additional complexity of including such a layer would obfuscate the processes examined in this study. Here, we strive to isolate the governing processes of instability development within the plume and surrounding mantle. Future work should examine the role of a (growing) thermal boundary layer on UMP.

4.5 Numerical resolution and upper mantle pulsing

Simulating long-lived upper mantle pulsing requires a well-resolved plume in the upper mantle. To test the resolution dependence, we ran 4 models with parameters similar to those presented previously (composite rheology, 410-km Clapeyron slope of 4 MPa K⁻¹, overlying plate speed of 8 cm yr⁻¹, and initial plume temperature anomaly of 550 K). We then set the minimum upper mantle element sizes to either 32, 16, 8 (reference resolution) or 2 km. In the poorest resolved case (32 km wide elements), no TZP occurs and UMP initiates at ∼360 Myr MT. For the same parameters, more finely resolved models (element size <32 km) successfully simulate both TZP and UMP (Type 1). All resolutions with element sizes <32 km exhibited similar UMP and/or TZP initiation times, plume conduit tilts and periodicities; TZP initiated in all these runs between 165 and 172 Myr, and UMP between 220 and 226 Myr. This suggests that element sizes of ≥32 km numerically inhibit plume instabilities, even when reaching large conduit tilts. Thus, to accurately resolve TZP and UMP models require elements <32 km across with 16-km-wide elements sufficient to capture the timing and periods of UMP in our models.

4.6 Comparisons to previous work

Although early analogue models demonstrated diapir formation from instabilities in a tilted plume conduit (Whitehead 1982), later studies suggest that such instabilities will fail to form in thermal plumes under mantle conditions (Kerr et al. 2008; Mériaux et al. 2011). The primary arguments against instability formation in the mantle include: (1) thermal diffusion helps stabilize low Rayleigh number plumes (Kerr et al. 2008) and (2) the theoretical growth rate of instabilities, in most cases, is larger than the time for material to rise through the upper mantle (Mériaux et al. 2011). However, these previous studies examine the rise of horizontal cylinders of buoyant, less viscous material in an initially motionless ‘mantle’ with either constant or temperature dependent rheology. No previous study has examined instability growth of tilted, buoyantly rising plume conduits in a mantle described by a realistic rheology, including phase changes, and with tilting driven by motion across the top of the domain.

Despite the above arguments against instability formation in tilted plume conduits at mantle conditions, our simulations exhibit instability growth for a range of plate speeds, plume temperatures, and rheological parameters (Section 3). To better understand instability formation in our models, we quantify the instability growth rate relative to plume rise rate by assuming that tilted plume conduits act similarly to a Rayleigh–Taylor instability. As a first order estimate of instability formation, we use an equation for instability growth rate appropriate for two layers of constant, but not identical, viscosity (eq. 10 from Houseman & Molnar 1997). The calculated instability growth rate is a function of plume radius, pulse wave number (⁠|${{\bf \lambda }}\ = \ 2\pi /( {P*v} )$|⁠), where v is the plate speed and P is the measured pulse period), and viscosity ratio between the plume and ambient mantle. We find that the resulting characteristic instability growth time (inverted value of the calculated growth rate) is less than the time to rise across the upper mantle through a tilted circular pipe (e.g. Turcotte & Schubert 2013; with tilt accounted for by multiplying the vertical pressure gradient driving flow by the cosine of the plume tilt) for viscosity ratios below a critical value, implying that instabilities should form for these viscosity ratios (Fig. 9a). Conversely, our calculation predicts that above a critical viscosity ratio the growth time will exceed the rise time and instabilities will fail to form in the upper mantle. The critical viscosity ratio depends on the plume conduit radius, the plume tilt, and the pulse wave number, and thus varies depending upon the plume parameters (Fig. 9b). We note that using the stokes velocity of a sphere (e.g. Mériaux et al. 2011) or that of a tilted conduit (Happel & Brenner 1965; Skilbeck & Whitehead 1978) are similarly appropriate for estimating rise times, however, neither of these methods accounts for both the viscosity ratio and plume tilt that our results suggest are important for instability development. Additionally, we find that the velocity calculated using the method presented here exhibits the fastest rise times and the smallest range of viscosity ratios for instability development (e.g. Fig. S5), providing the most restrictive estimates of the critical viscosity ratio for instability development. Based upon this analysis, generally lower wave numbers (longer periods or higher plate speeds) and smaller plume radii promote instability growth by increasing the instability growth rate (decreasing the growth time, Fig. 9b) and increasing the rise time. Plume tilt also affects the likelihood of plume instability development. Increasing the plume tilt increases the rise time making instability development more likely. For example, given a plume radius of ∼50 km and wave number of ∼4 × 10⁻⁵ m⁻¹ increasing the tilt from 33.5° to 67° increases the critical VR from ∼20 to ∼60, approximately doubling the range of VRs predicted to develop instabilities and start UMP. Note that at very low wavenumbers (<∼1 × 10⁻⁵ m⁻¹) the likelihood of instability development rapidly decreases regardless of plume radii or plume tilt.

Figure 9.

(a) Rise and instability times relative to the viscosity ratio (VR, calculated from eq. 10 in Houseman & Molnar 1997) between the plume and ambient mantle for fixed values of plume radius and wavenumber (as indicated). The black dot represents the critical VR where instability time first exceeds rise time indicating the maximum VR where instability formation is predicted to be possible. In this example, the VR must be below ∼40 for instabilities to form. (b) The predicted critical VR (contours) value varies considerably for different plume radii and wavenumbers, and plume tilt (solid black lines: 67°; dashed cyan lines: 33.5°; solid, violet lines: 0°). The grey box encloses the approximate range of radius/wave numbers present in the plumes modelled in this study. The red dot indicates the reference model's approximate plume radius and instability wavenumber at the time of the earliest periodicity taken using the running window (e.g. Fig. 5), and the tip of the red arrow indicates the final periodicity from the running window to showcase the trend of plume radius and wavenumber through model evolution.

Open in new tabDownload slide

Mériaux et al. 2011 suggest that 3-D plume conduits should be more stable and less prone to tilting related instabilities than 2-D plumes. Although a full 3-D study is beyond the scope of this paper, we perform a restricted test of plume pulsing by simulating a rising mantle plume in a 3-D box with composite rheology and dimensions of 2000 × 400 × 600 km (X, Y and depth), representing a section of the upper mantle and transition zone. The plume is initiated by prescribing a Gaussian function of temperature and velocity on the bottom boundary (i.e. injecting warm fluid at the base of the domain), with a maximum amplitude of 25 cm yr⁻¹ and a purposefully rapid overlying plate speed of 20 cm yr⁻¹ (for other model parameters, see the parameter, i.e. prm, file in zenodo). In this simplified example, the plume is more stable than the 2-D cases, however, upon reaching a tilt of ∼90° pulses begin to form in the plume stem. Although very preliminary in nature, this simulation suggests that 3-D plume pulsing may require steeper tilts relative to 2-D plumes, but that pulsing may still be possible at mantle conditions. However, more extensive models are needed to fully test pulsing of tilted 3-D plume conduits.

4.7 Comparisons to real-world plumes

4.7.1 Timescale for initiation of plume pulsing

UMP initiates >51 Myr (Table 4) after a plume head arrives in the upper mantle (1–21 Myr after reaching maximum tilt), suggesting that only the longest hotspot volcanic chains currently visible on our planet could record this type of igneous variability. Such hotspot chains include Hawaii, Louisville, and Kerguelen (Ninety East Ridge), of which observations at both Hawaii and Louisville show periodic melt production at periods similar to our models (∼2–8 Myr; Van Ark & Lin 2004; Vidal & Bonneville 2004; Wessel 2016; Morrow & Mittelstaedt 2021). However, in our simulations plume tilt is caused solely by an imposed surface motion and the resulting shear flow in the upper mantle. We expressly avoid other processes likely to cause plume tilt to isolate the evolution of individual plumes. Yet, 3-D spherical numerical models of the entire mantle show plumes drifting at relatively rapid rates of 1–5 cm yr⁻¹ due to additional processes including slab push and plume merging (Arnould et al. 2020). If such drift causes plume tilting in the upper mantle, a plume could reach our proposed critical tilt angle of 67° in 10–50 Myr and, thus, UMP initiation as quickly as 11 Myr after a plume reaches the upper mantle. Combining the shear flow in our models with plume drift due to the above processes could cause even more rapid plume tilt, suggesting that hotspot tracks with much smaller age extents may also record the effects of UMP if it occurs.

4.7.2 Plume periodicity and viscosity ratio

In our numerical simulations, we observe pulsations of tilted plume conduits resulting in variations in plume volume flux over periods of (primarily) 1.6–6.7 Myr (Table S2). Although studies examining variability in plume melt production (a proxy for mantle flow variability) suggest that the upwelling rate of individual plumes may exhibit periodic variations over a wide range of timescales (e.g. Van Ark & Lin 2004; Mjelde & Faleide 2009; Sreejith & Krishna 2015), a recent study of twelve hotspot volcanic chains finds eight hotspots (∼67 per cent) with at least one statistically significant period of melt production variability in the ∼1.6–6.7 Myr range (Morrow & Mittelstaedt 2021). Data from hotspots not included in the above study, such as Amsterdam St Paul (∼6 Myr; Maia et al. 2011; Sibrant et al. 2019) and Iceland (∼8 Myr, Parnell-Turner et al. 2014), also show evidence of similar periods of melt flux variability. Additionally, Adam et al., 2022 suggest that ˜5-10 Myr signals observed from estimates of magmatic and swell flux variations at the Louisville hotspot are likely caused by shallow phenomena, in particular to instabilities within tilted plume conduits.

For UMP to occur, our results and scaling analysis (Fig. 9) suggest that the VR between the ambient mantle and a tilting plume must be less than a critical value. Thus, if UMP is responsible for observed igneous variability, we expect real-world plumes to exhibit a similar VR limit. Taking the plume excess temperatures estimated by Morrow & Mittelstaedt (2021), assuming an ambient upper mantle potential temperature between 1585 K (1312 °C; Ball et al. 2021) and 1727 K (1454 °C; Putirka et al. 2007), a strain rate within the plume stem of 1 × 10⁻¹³ s⁻¹, and the same rheology as our models (eqs 16–18), we estimate that all 12 hotspots examined by Morrow & Mittelstaedt (2021) have a VR <∼20. This value is well within our predicted VR limit for UMP at all but the largest wave numbers and plume conduit radii tested. However, if we use a lower strain rate estimate of 1 × 10⁻¹⁴ s⁻¹ and the cooler estimate for mantle potential temperature (1585 K), only plumes with excess temperatures ≤205 K will have a VR < 20. VR values for warmer plumes such as Hawaii and Galápagos (ΔT ≈ 230 K; VR ∼ 35) would predict UMP to require plume conduit radii ≤50–60 km. Given the combined uncertainties on plume excess temperature (e.g. table A1 in Morrow & Mittelstaedt 2021), mantle potential temperature (e.g. Putirka 2005; Putirka et al. 2007; Katsura et al. 2010; Dalton et al. 2014; Ball et al. 2021), and mantle viscosity, it is difficult to tightly constrain the VR of real-world plumes. Even so, the above calculations suggest many plumes may have VRs consistent with those predicted to be required for UMP.

4.7.3 Plume tilt

In addition to measurements that yield estimates of pulsing periodicity and mantle to plume VR, the tilt of plume stems could provide a useful comparison to our models. Indeed, interpretations of seismic tomography suggest that several plumes are tilted including Hawaii (French & Romanowicz 2015), Galápagos (Villagómez et al. 2014), St Helena (French & Romanowicz 2015), Tristan (Schlömer et al. 2017), Iceland (Celli et al. 2021) and Reunion and Bouvet (Tsekhmistrenko et al. 2021). Despite being able to infer that a plume is tilted in current tomographic models, quantifying the exact tilt of these plumes is still difficult at current tomographic resolutions.

4.7.4 Pulse melt production

Changes in plume volume flux in the uppermost mantle will cause approximately synchronous changes in hotspot melt production. Numerically, the average deviation in 2-D melt flux by a plume pulse varies from 3.8 × 10⁻⁵ to 26 × 10⁻⁵ km³ yr⁻¹ km⁻¹. Using these average melt fluxes, we can calculate a thickness of material produced (per km perpendicular to our model) as the plate moves overhead at 8 cm yr⁻¹. Assuming a pulse duration is 3–5 Myr (indicating the plate moves 240–400 km), we predict an upper mantle pulse to produce excess melt equivalent to a constant layer of material ∼2 km thick (∼0.48–3.2 km).

For comparison, we can use the above estimate in two simplified end-member scenarios: (1) all the excess melt is erupted as a layer with uniform thickness during a pulse and (2) all the excess melt is emplaced as a uniform layer of underplating at the base of the crust. In the case where excess melt erupts to the seafloor, adding 2 km to the thickness of a hotspot volcano should be readily visible in bathymetry; this is equivalent to ∼20 per cent of the seafloor to summit height of Mauna Loa volcano on the Big Island of Hawaii. We note that there are many caveats to this comparison (e.g. elastic plate flexure, melt transport, eruption from a centralized vent, etc.) and only provide this comparison to place our calculated melt volumes in a relatable context.

Alternatively, if we assume that a pulse will emplace an additional ∼2 km of underplating beneath the hotspot track, we can compare this thickness to seismic and gravity constraints on underplating and the isostatic contribution such a layer would make to the hotspot swell. Underplating at Hawaii is suggested to be ∼4–7 km thick based upon seismic receiver functions (Leahy et al. 2010) or up to 10 km thick based upon gravity analyses (Morrow & Mittelstaedt 2021). Thus, a 2-km-thick layer of underplating is equivalent of ∼20–50 per cent of underplating thickness estimates beneath the Hawaiian chain. Next, considering the isostatic response of the crust, a 2-km-thick layer with a 300 kg m⁻³ density deficit relative to the mantle will add ∼250 m to the hotspot swell, 25 per cent of the observed along-track swell variability at Hawaii (∼1 km; Adam et al. 2005). If we assume a 3-D plume to have Poiseuille flow one would expect a higher volume flux and melt production in the conduit centre, which may lead to higher intrusive and/or extrusive thicknesses near the hotspot centre, increasing the expected swell and/or volcano volume changes. However, it is unclear if this assumption is valid since recent work accounting for grain size evolution suggests that plug flow (uniform upwelling velocity across the conduit) may be more likely within a plume stem (Dannberg et al. 2017). Overall, we suggest that our models predict that UMP would produce measurable changes in hotspot melt production, regardless of whether that melt was erupted or intruded.

4.8 Model limitations and future work

This study examines the dynamic evolution of a whole-mantle, 2-D thermal plume with realistic rheological behaviour and major mantle transition zone structure subject to a constant lateral plate velocity at the model surface. Results suggest that UMP could explain one facet of variability in plume flux over time. However, observations find many hotspots show a larger range of periodicities (e.g. Morrow & Mittelstaedt 2021), suggesting that if UMP occurs, it is unlikely to be the only process influencing plume dynamics and, thus, volcanic flux. We neglect several processes in our models that may alter the periodicity of UMP such as interaction with an overlying TBL, grain size evolution or compositional heterogeneity. Additionally, consideration of the third dimension could lead to increased stability of plumes (Mériaux et al. 2011), perhaps impeding UMP development.

Grain size evolution may play a large role in the viscosity structure of upwelling plumes and, thus, plume dynamics. Unlike simulations with diffusion creep and constant grain size where the lowest viscosities are located at the centre of the plume stem, simulations with dynamic grain size predict lower stresses in the plume centre to cause increasing grain sizes, raising viscosity (Dannberg et al. 2017). This is compounded by high stresses on the edges of the plume that reduce grain size and decrease viscosity (Dannberg et al. 2017). Overall, grain size evolution within plume stems predicts viscosities like those predicted by dislocation creep, where high strain rates result in reduced viscosity (Dannberg et al. 2017). Unlike dislocation creep, however, the effects of grain size evolution depend strongly on the rate of change of strain rate instead of its magnitude (Dannberg et al. 2017). Importantly for UMP development, evolving grain sizes can lead to increased decoupling of the upper and lower mantle, and consequently greater plume tilts (Dannberg et al. 2017).

In addition to grain size evolution, compositional variability within a plume conduit, perhaps due to entrained material, can alter plume dynamics (Kumagai et al. 2008). This study considers compositionally homogenous thermal plumes, but isotopic data from plumes such as Hawaii suggest some plumes may vary compositionally and spatially over time (e.g. Weis et al. 2011). Furthermore, studies indicate that compositional differences such as the entrainment of eclogite into a plume can lead to pulsations of plume material in the upper mantle at rates comparable to some of the variability in Hawaii (Ballmer et al. 2013).

5 CONCLUSIONS

This paper uses numerical models to examine the formation of shallow instabilities, or pulses, in tilted plume conduits as a potential mechanism for changing hotspot melt production through time. Our 2-D models simulate 400 Myr of the evolution of a single, whole-mantle plume as it rises through the mantle and spreads beneath a moving plate. To simulate plume dynamics, our models solve the conservation equations for a compressible mantle subject to a temperature-, and depth-dependent linear rheology or a temperature-, depth-, and strain-rate dependent non-linear rheology that includes a thermodynamically consistent treatment of phase changes in the mantle transition zone. We evaluate model results by quantifying changes in upper mantle vertical plume flux, plume tilt and melt production for different values of Clapeyron slope, horizontal surface velocity (plate motion) and plume temperature.

Our model results indicate that plume conduits that tilt >67° from vertical develop instabilities that lead to repeated plume pulsing in the upper mantle [i.e. upper mantle pulsing (UMP)]. UMP initiates after the plume tilts in one of two ways: Type (1) plate shear drives and maintains the plume tilt and Type (2) deeper pulsing from the transition zone continues into the upper mantle, destabilizing the tilted plume stem. In all models where UMP occurs, after initiation pulsing continues until the end of a model run (400 Myr).

Frequency analysis of plume volume flux in our models demonstrates that UMP forms regular pulses with statistically significant (99 per cent confidence) periods. The number and distribution of pulsing periods depend on the model rheology, where more effective non-linear rheology results in fewer significant signals and an overall reduction in pulse periodicity. Regardless of rheology, models primarily yield a few closely clustered periods, with a range from ∼2 to 6 Myr (Fig. S4; Table S2). Results indicate that pulsing period depends on the viscosity ratio (VR) between the plume and ambient upper mantle, with lower VRs promoting pulsing by decreasing the instability growth time relative to plume rise time through the upper mantle. Additionally, analysis of plume volume flux after initiation of UMP using a 25 Myr wide moving window reveals that pulsing period increases as VR decreases and overall mantle viscosity increases.

Based upon comparison with observations, we conclude that UMP may contribute to plume melt production variability along hotspot volcano chains. Although UMP in our models initiates >51 Myr after a plume arrives in the upper mantle, UMP type 2 suggests that other processes (such as TZP in the models presented here) can destabilize a plume and initiate UMP. As such, additional mechanisms for tilting plumes could reduce this time to well within the duration of many hotspot island chains (>∼11 Myr). Furthermore, estimates of changing melt production at several hotspots find periods of ∼5 Myr, consistent with the range of periods observed in our models (∼1.6–6.7 Myr).

Further work is needed to predict whether UMP is a feasible explanation for variations in magmatic output along hotspot tracks. Such work should include fully 3-D numerical simulations with appropriate boundary conditions and sufficient resolution (<32 km element widths) to resolve UMP (Section 4.5). Additionally, the impact of grain size evolution on the rheology of upwelling plumes may alter formation of UMP.

SUPPORTING INFORMATION

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials sup- plied by the authors. Any queries (other than missing mate- rial) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This study was supported by NSF awards EAR-1520856 and OCE-1753354 to EM. We thank the Computational Infrastructure for Geodynamics (geodynamics.org), which is funded by the National Science Foundation under award EAR-0949446 and EAR-1550901, for supporting the development of ASPECT. Models were run using the Computational Resources Core of the Institute for Interdisciplinary Data Sciences (IIDS) at the University of Idaho, or using the ETH Zurich Euler cluster. Figures in this paper were made using Paraview, InkScape and Matlab and colourscales from Crameri (2018). We would like to thank John Naliboff for his help in initially setting up the ASPECT models, and John Naliboff and Jerry Fairley for their comments on the thesis this study builds upon. We would also like to thank Bernhard Steinberger and an anonymous reviewer for their very detailed and constructive reviews.

DATA AVAILABILITY

Models for this study were conducted using Dealii version 8.5.1, and ASPECT version 2.0.0-pre using a modified material model and adiabatic conditions to include Clapeyron slope related temperature jumps. The specific ASPECT branch and parameter files to reproduce these models, as well as the scripts used to create figure 9 and calculate the volume flux, periodicity, plume tilt, melt, and volume flux can be found at https://zenodo.org/record/7299391. The author contributions to this manuscript are as follows: DN ran the models, analysed the data, interpreted results, and cowrote the manuscript, and EM conceptualized the project, helped with model set-up, aided in the interpretation of results and cowrote the manuscript.

REFERENCES