A joint geochemical–geophysical record of time-dependent mantle convection south of Iceland

The North Atlantic V-Shaped Ridges (VSRs) provide a spatially extensive and clear record of unsteady mantle convective circulation over > 40 My. VSRs are diachronous ridges of thick crust formed with a periodicity of ∼ 5 My along the Mid Atlantic Ridge, south of Iceland. We present data from a set of dredged basalt samples that shows chemical variation associated with two complete VSR crustal thickness cycles where they intersect the Mid Atlantic Ridge. The new dataset also records chemical variation associated with a VSR crustal thickness cycle along a plate spreading ﬂow-line. Inverse correlations between crustal thickness and both incompatible trace element concentrations and incompatible element ratios such as Nb/Y and La/Sm are observed. Geochemical and crustal thickness observations can be matched using a time-dependent mid-ocean ridge melting model with a basal boundary condition of sinusoidally varying potential temperature. Our observations and models suggest that VSRs are generated when hot patches are carried up the plume stem beneath SE Iceland and spread radially outward within the asthenosphere. These patches are then drawn upward into the melting region when passing beneath the Mid Atlantic Ridge. The geometry of the VSRs and the size of the dynamically supported swell suggest that the Iceland Plume is the strongest plume in the Earth at present, with a volume ﬂux of 49 ± 14 km 3 yr − 1 . © 2013 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).


Introduction
Convection within the upper mantle is expected to be timedependent because the Rayleigh number is super-critical by 3 to 5 orders of magnitude (Schubert et al., 2001). Time-dependence can take two basic forms, one in which convection cells move laterally with respect to one another (plume wander) and another in which patches of differing temperature are advected round a convection cell (plume pulsing) (White and McKenzie, 1995). The North Atlantic V-Shaped Ridges (VSRs) provide a long period, spatially extensive and clear record of the plume pulsing type of unsteady mantle convection over time periods of order 1 to 10 million yr (Vogt, 1971;Jones et al., 2002b). VSRs are diachronous ridges of thick crust formed at the Mid Atlantic Ridge to the north and south of Iceland (Fig. 1). Since their discovery, it has been generally agreed that the diachronous geometry results from melting anomalies that propagate away from Iceland within the asthenosphere (Vogt, 1971;Jones et al., 2002b;Poore et al., 2011). The process that generates VSRs also appears to modulate Atlantic oceanic circulation, since the VSR record correlates with stable isotope proxies for meridional overturning circulation (Wright and Miller, 1996;Poore et al., 2006).
The North Atlantic VSRs comprise one of the most important mantle convection records on Earth because no other plume-ridge system records as many melt pulses over such a large distance from a hotspot. However, knowledge of geochemical variability associated with the VSRs is lacking in comparison with available geophysical records of crustal thickness and oceanographic records of deep-water flow. There is still debate over whether VSR melt fluctuations are caused by thermal or compositional variability in the mantle, or whether they reflect neither and result instead from crustal accretion processes (Jones et al., 2002b;Hey et al., 2010). Here we report dredged basalt samples obtained during RV Celtic Explorer cruise CE0806.  (a) Gravity proxy for oceanic crustal structure: high-pass filtered free-air gravity field, after Jones et al. (2002b). Dredge locations with black centres yielded samples of local basalt. Plume centre from Shorttle et al. (2010). Coloured arcs show traces of hot and cool rings within the asthenosphere, and correspond to vertical coloured stripes in Figs. 3, 4 and 6. (b) Relationship between 2D model slice and 3D ridge-plume interaction geometry, showing likely path of mantle arriving and melting beneath study area. a record of basalt geochemistry covering two complete VSR crustal thickness cycles where they intersect the Mid Atlantic Ridge, as well as a VSR crustal thickness cycle along a plate spreading flowline.
Coincident geochemical and geophysical datasets are key to determining the cause of VSR melting anomalies. We analyse seismic and gravity data to estimate crustal thicknesses along V-shaped ridge crests, V-shaped troughs and the Mid-Atlantic Ridge axis (Section 3). We identify the geochemical signature of the VSRs by correlating geochemical proxies for degree of melting with crustal thickness (Section 5). We then develop a numerical melting model that tracks passage of hotter and cooler mantle patches through the melting region beneath a mid-ocean ridge axis (Section 6). This model can reproduce the main characteristics of both geochemical and crustal thickness VSR datasets. The new model provides estimates of the temperature variation between hotter and cooler mantle patches, rather than the average temperature change within the melting region estimated by previous studies Poore et al., 2011). Finally, we discuss implications of our results for mantle plume flux (Section 7).

V-Shaped Ridges and the Iceland Plume head
VSRs are topographic seafloor ridges built of mid-ocean ridge basalt that were originally recognized on single-channel seismic profiles (Vogt, 1971). Here we exploit the satellite-derived gravity field as a proxy for VSR topography (Fig. 1a). This proxy is helpful because it provides an extensive map of uniform quality and also because it sees through sedimentary cover up to 1 km in thickness (Jones et al., 2002b). Red VSR stripes on Fig. 1a represent topographic ridges supported by thicker crust and blue VSR stripes represent troughs underlain by thinner crust. Gravity, topography and crustal thickness are positively correlated in this region (Navin et al., 1998;Smallwood and White, 1998;Poore et al., 2009). Another set of stripes caused by fracture zones can be seen on Fig. 1a. The gravity and topography signatures of these fracture zones have higher amplitudes than the corresponding VSR signatures. Therefore gravity and topography maps cannot be used to determine whether VSR crustal thickness anomalies exist in the fractured crust (Jones et al., 2002b). Vogt's suggestion (1971) that nested ridges of thicker and thinner crust are produced when melting anomalies propagate outward from Iceland is widely accepted, but the cause of the propagating melting anomalies is still debated. One hypothesis is that plume head asthenosphere travelling outward from Iceland contains patches of differing temperature that are drawn upward into the melting region beneath the Mid Atlantic Ridge Jones et al., 2002b;Murton et al., 2002;Poore et al., 2011). Alternatively, the patches within the plume head could differ in lithology (Foulger et al., 2005). In these two hypotheses, thermal or compositional patches rise up the plume conduit before spreading within the plume head. A third possibility is that flow up the plume conduit is steady and fluctuations in flux within the plume head occur when the flow at the top of the conduit interacts with a stepped lithospheric base inherited from rift relocations (Harðarson et al., 1997;Sleep, 1996). Most recently, Hey et al. (2010) and Benediktsdóttir et al. (2012) have suggested that the VSRs do not result from melting anomalies, but instead represent a series of rift propagation cycles. A joint geochemicalgeophysical record is required to test these hypotheses.
VSR melting anomalies travel within the head of the Iceland mantle plume. We have argued that plume head flow occurs radially outward within the asthenosphere and is not significantly channelled beneath the spreading axis (White and Lovell, 1997;Jones et al., 2002b;Shorttle et al., 2010;Poore et al., 2011). Others have suggested that plume-driven flow is preferentially channelled beneath the Mid Atlantic Ridge (Vogt, 1971;Albers and Christensen, 2001;Gaherty, 2001;Tilmann and Dahm, 2008). Nevertheless, observations that support radial flow are compelling. First, the Iceland plume swell is roughly circular (Fig. 2). Secondly, this dynamically supported region spans several thousand kilometres from continental Greenland to the continental margins of NW Europe, whereas the melting region beneath the spreading ridge is <100 km wide (Section 6). Thirdly, a seismically slow patch within the asthenosphere is imaged beneath the entire dynamically supported swell, including the edges beneath the adjacent continents (Pilidou et al., 2005;Rickers et al., 2013). Fourthly, joint modelling of trace element compositions and crustal thickness observations in central Iceland indicates that little plume-driven flow occurs shallower than about 100 km (Maclennan et al., 2001). Finally, correlation between Miocene-recent deep-water flow and VSR records indicates that plume-head processes affect vertical motion of the Greenland-Scotland Ridge over 500 km away from the spreading axis (Wright and Miller, 1996;Poore et al., 2006Poore et al., , 2011Robinson et al., 2011).
The plume centre lies in southeastern Iceland (Shorttle et al., 2010). The position of the plume centre has been determined from the radial symmetries of the plume swell and crustal thickness anomalies, as well as regional and global seismic studies. The depth from which the plume rises is debated (Foulger et al., 2001;Rickers et al., 2013) but this issue does not affect the results of this study. The geochemical and geophysical observations we discuss depend on the temperature structure that results from interaction between mid-ocean ridge spreading and the top of the plume head far from the conduit. For our purposes, the plume conduit can be regarded as a point source that lies within the asthenosphere at the centre of the plume head.

Crustal thickness signature
Both the thickness and composition of oceanic crust are required to test detailed models of VSR melting. In particular, we require estimates of the difference in crustal thickness between V-shaped ridges and adjacent troughs. Three wide-angle seismic experiments constrain crustal thickness in our study area (Fig. 3a). Previous studies interpolated a single line through these data to predict crustal thickness as a function of distance from Iceland. Here, we divide the crustal thickness measurements into two groups, representing V-shaped ridge crests and V-shaped troughs (Fig. 3b). The ridge group contains three measurements: a point measurement near the crest of VSR-2 (Bunch and Kennett, 1980); VSR-2 near the eastern end of the RAMESSES profile (Navin et al., 1998); and VSR-1 seen on profile CAM-74 (Smallwood and White, 1998). The trough group contains two measurements: the intersection of profile RAMESSES with the MAR close to VST-1 (Navin et al., 1998); and the intersection of profile CAM-74 with VST-1 (Smallwood and White, 1998). These seismic experiments provide reliable constraints on crustal thickness because they recorded Moho reflections, turning rays within the lower crust and (except Bunch and Kennett, 1980) turning rays within the mantle. All of the seismic studies carried out resolution tests which suggest that the crustal thickness measurements are accurate to between ±0.5 and ±0.8 km.
Separate model lines were fitted to the ridge and trough data groups. Model lines were constrained at their northern and southern ends by wide-angle seismic profiles immediately south of Iceland and near the Charlie Gibbs fracture zone (Weir et al., 2001;Whitmarsh and Calvert, 1986). Crustal thickness along the Mid Atlantic Ridge is expected to vary between the ridge and trough endmember curves. Detail of this variation is provided by the gravity proxy for local crustal structure (Fig. 1). The bottom panel of Fig. 4 shows the mean of the gravity proxy field on either side of the spreading axis (Jones et al., 2002b). Averaging the conjugate gravity fields reduces noise associated with asymmetric crustal accretion (Hey et al., 2010;Benediktsdóttir et al., 2012) and enhances the VSR gravity signature thought to be generated by mantle melting. The mean field was projected parallel to the VSR crests to obtain the gravity curve in Fig. 4. This projection substantially reduces noise associated with fracture zones and further reduces noise associated with asymmetric accretion. A crustal thickness model for the spreading axis was then obtained by stretching the projected gravity curve between the ridge and trough end-member curves (Fig. 3).

New geochemical dataset
RV Celtic Explorer cruise CE0806 (April-May 2008) carried out 47 dredges at 33 locations between 55 • and 62 • N along the axis and western flank of the Mid Atlantic Ridge (MAR) over a period of 12 days (Fig. 1a). 37 dredges returned local basaltic rock suitable for chemical analysis. The remainder consisted of coral and glacial dropstones of major Greenland and minor Icelandic provenance. Basaltic dropstones were recognized by their smooth surfaces and  Jones et al. (2002a). Residual topography obtained by correcting GEBCO ocean depth compilation for sediment loading using sediment thickness compilation of Louden et al. (2004) and subtracting age-depth curve (100 km thick plate) of Crosby et al. (2006). Circle X, plume head radius of ∼500 km assumed by several studies when calculating Iceland Plume flux (Section 7); circle Y, minimum plume head radius of 1370 km from new geochemical dataset; circle Z, plume head radius of 1800 km estimated from dynamic topography. Arrows indicate Greenland Scotland Ridge. (b) Trace element and isotopic compositions measured along the Mid Atlantic Ridge within the blue circle. Additional data from Schilling et al. (1983Schilling et al. ( , 1999, Devey et al. (1994), Murton et al. (2002). lack of glassy rind and excluded from chemical analysis. Dredge sites are grouped into 3 transects. Transect 1 collected very young rock along the MAR axis between 55 • and 57.5 • N. The aim was to measure chemical variation associated with the intersection of the MAR and the crest of VSR-2 and the adjacent trough VST-2. This dataset augments a set of samples collected along the MAR between 57.5 • N and the southern edge of the Iceland Shelf (Murton et al., 2002). The combined dataset crosses the intersections of the MAR with VSR-1, VST-1, VSR-2 and VST-2 (i.e. 2 complete VSR ridge-trough cycles). Transects 2 and 3 sampled the western arms of VSRs 2 and 1 respectively. These transects were designed to compare the compositions of V-shaped ridges and troughs formed along plate spreading flow-lines.
XRF analyses on whole-rock samples and glass chips were carried out at the University of Edinburgh using techniques similar to those described by Fitton et al. (1998) with modifications noted by Fitton and Godard (2004). ICP-MS analyses and isotopic analyses were carried out at the National Oceanography Centre, Southampton using procedures similar to those described by Murton et al. (2002). Further details of analytical procedures are provided in Section A (Supplementary Material). Accuracy was determined by analysis of USGS and GSJ geological standards. Precision was estimated by repeat analysis of international standards and a local sample from the Reykjanes Ridge. A glassy basalt from the south-ern Reykjanes Ridge (BRR-1) was re-analysed to ensure consistency between the Murton et al. (2002) dataset and the new dataset. A complete set of chemical analyses, including standards, together with accuracy and precision estimates is provided as a Supplementary Dataset.

Geochemical signature
Geochemical data and crustal thickness estimates along the Mid Atlantic Ridge (MAR) between 500 and 1400 km from the plume centre (56 • to 61.5 • N) are plotted in Fig. 4. The study area lies south of the region strongly enriched in incompatible trace elements and radiogenic isotopes surrounding Iceland (Fig. 2b). We seek relatively small geochemical signals that correlate with VSR crustal thickness variations.
Incompatible element concentrations (e.g. La, Nb) and incompatible element ratios (e.g. La/Sm, Nb/Y) show several components of variation with distance from Iceland. An almost-periodic component on the ∼400 km scale can be seen clearly on uninterpreted data. This component is superimposed upon longer wavelength variation. In order to quantify the periodic component of the variation we used both the non-parametric regression technique of Samworth and Poore (2005) and also the best-fitting order 6 polynomial curve. Order 6 was chosen because four MAR/VSR intersection are expected and there are two end points. The polynomial  Bunch and Kennett (1980); Nea98, Navin et al. (1998); SW98, Smallwood and White (1998). Measurements assigned to ridge (red) or trough (blue) and interpolated using akima splines to give end-member crustal thickness models. Thick black line: new crustal thickness model for Mid Atlantic Ridge axis with uncertainty bound from ±1 standard deviation uncertainty on the projected gravity field; thin black line, Jones and Maclennan (2005); dashed line, . Vertical bands: VSR/axis intersections from Fig. 4. regression lies within the confidence band of the non-parametric regression, and corresponding extrema lie at similar distances from Iceland.
Gravity data show that two V-shaped ridges (VSR-1 and VSR-2) and two V-shaped troughs (VST-1 and VST-2) intersect the MAR in the study area. (VSR naming conventions in this and other papers are compared in Table S1, Supplementary Material). The MAR/VSR-1 and MAR/VST-1 intersections are seen clearly because the crust is not dissected by fracture zones. Crust further out from Iceland is fractured, so that the MAR/VSR-2 intersection is less clear and the MAR/VST-2 intersection cannot be directly seen. VSR crustal thickness variations could be present in the fractured crust (Section 2). Therefore, models for the geometry of all four VSRs under consideration have been fitted to the gravity data to clarify the points of intersection between the MAR, VSR-2 and VST-2. Model parameters are listed in Table S2 (Supplementary Material) and the modelling procedure is described in Section D (Supplementary Material).
Coloured vertical stripes on the right-hand panel of Fig. 4 are centred on MAR/VSR intersection points. These stripes highlight correlation between the periodic components of compositional variation and crustal thickness. Lower trace element concentrations and ratios correspond to higher crustal thickness. For a given VSR, peak or trough positions estimated from different data series are spread over several 10s of kilometres. Widths of the vertical stripes marking the VSR/MAR intersection points were chosen by eye to encompass the lateral spread in extrema.
Trace element variation associated with VSRs cannot be explained by fractional crystallization because the magnesium number (Mg#) shows no corresponding periodic component of variation as a function of distance. Correlation between Mg# and residual scatter about non-parametric regression lines on Fig. 4 is observed for all trace elements but not for trace element ratios (Fig. 5). This correlation results from fractional crystallization. Murton et al. (2002) estimated that the extent of crystallization is 10-15%. The fractional crystallization trends for sample subsets with Nb/La < 0.8 and 0.8 < Nb/La < 1.2 are similar to trends obtained for the entire sample set, and there are few samples with larger Nb/La. There is therefore no evidence in this dataset that melts of different initial compositions have significantly different fractional crystallization behaviours. Fractional crystallization accounts for 15-40% of the total variance in scatter about the non-parametric regression lines on Fig. 4. Analytical uncertainty accounts for only a few percent of the total variance. We interpret that the remaining 60-80% of the variance in scatter comes from incomplete mixing of melts of different initial compositions. There is no evidence that these melt mixing or emplacement signals are correlated with the periodic VSR signal observed on Fig. 4. We conclude that the VSR-related trace element variation is a signal from the mantle.
The Sr isotope curve shows fluctuations that correlate fairly well with crustal thickness. No fluctuation in εNd corresponding to VSR-1, VST-1 and VSR-2 is discernible. There is abundant evidence for small-scale heterogeneity in the mantle source from various isotopic systems and trace element-isotope correlations (Murton et al., 2002;Kokfelt et al., 2006;Maclennan, 2008;Shorttle et al., 2010). Fluctuations in Nd and Sr isotopes that correlate with the VSR-related trace element fluctuations are therefore to be expected. Presumably the amplitude of any VSR-related isotopic fluctuations in εNd is smaller than that of the trace element indicators in Fig. 4 and is swamped by noise. We carried out principal component analysis on REE concentrations (Fig. 6). Only the first two principal components can be interpreted in terms of geochemical processes, and together they account for 97.5% of the data variance. The first principal component, PC 1 , accounts for 79% of the data variance. It has uniform weighting across the REEs, indicating a process that enriches or depletes all elements uniformly. A plot of PC 1 against radial distance resembles plots of Mg# or Yb concentration. We interpret that PC 1 reflects a combination of variation in degree of melting and fractional crystallization. The second principal component, PC 2 , accounts for 18.5% of data variance and has a negative slope from La to Lu, pivoting between Eu and Gd. This pattern is known to be associated with element fractionation by silicate melting, particularly when some melting occurs deep, in the presence of garnet (Slater et al., 2001;Maclennan et al., 2003;Fig. 4. Comparison of geochemical and geophysical VSR signatures. Bottom panels: gravity proxy for oceanic crustal structure. Gravity curve: mean ±1 standard deviation of gravity proxy projected along radial flow model lines. Crust curve from Fig. 3. Mg# is the Magnesium Number, defined as atomic Mg/(Mg + Fe 2+ ), and we assume that 90 mol% of the total Fe occurs as Fe 2+ . Error bars on geochemical data are 2σ ranges from repeat measurements of samples BRR1 (XRF) and BIR-1 (ICPMS). Right-hand panels: red lines, best-fitting degree 6 polynomial; light grey envelopes, confidence limits from non-parametric regression; vertical coloured stripes, interpreted VSR-Mid Atlantic Ridge intersections; details of model flow lines on gravity proxy field provided in Table S2 and Section D (Supplementary Material).  Rudge et al., 2013). A plot of PC 2 against radial distance shows a signal similar to the periodic component of variation in Nb, La, La/Sm and Nb/Y. We interpret that PC 2 corresponds to melting fluctuations associated with the VSRs. Principal component analysis therefore indicates that the VSR-related compositional fluc-tuations seen in Fig. 4 are a general feature of the whole trace element dataset.
The relationship between degree of melting and crustal thickness is an important indicator of the melting process (Fig. 7). Degree of melting was estimated using selected trace elements and ratios. Element concentrations were corrected for fractional crystallization using the regression lines in Fig. 5. A box-car filter of width 50 km was applied to simulate mixing of varying local instantaneous fractional melt compositions. This filter width corresponds to the typical distance between fracture zones. Crustal thickness was estimated for each sample using the model developed in Section 3. All filtered datasets show negative correlation between trace element indicator and crustal thickness, when viewed over the full crustal thickness range. This relationship indicates a positive correlation between degree of melting and crustal thickness. Some trace element indicators show a more variable melting-crustal thickness relationship at crustal thicknesses below 7.5 km. It is probably unwise to interpret this observation in terms of VSR formation processes, considering both scatter in the geochemical data arising from incomplete magma mixing and emplacement processes, and uncertainty in the crustal thickness model owing to the small number of wide-angle seismic constraints.
We conclude that we have isolated the geochemical signature of the North Atlantic VSRs. Geochemical and geophysical VSR signatures are strongly correlated. Both signatures likely reflect degree of melting, which suggests that the temperature of the mantle source is the principal control on the VSR chemical signature and on the origin of the VSRs themselves.

Modelling time-dependent melting
We have developed a model of time-dependent mid-ocean ridge melting to calculate volumes and compositions of erupted melts, in order to test the hypothesis that temperature of the mantle source is the principal control on VSR generation. This model uses standard calculation procedures and parameterisations to determine the thermal structure and melt productivity beneath a mid-ocean ridge, and it has two important new features. First, we apply a basal boundary condition of sinusoidally varying potential temperature to simulate the effect of successive pulses of hot mantle passing beneath the melting region within the plume head. Secondly, we track potential temperature and progress of melting throughout the melting region in order to calculate trace element compositions.

Model development
The model tracks thermal evolution of a two-dimensional section parallel to the direction of plate spreading. The relationship between the modelled region and the North Atlantic plume-ridge interaction scenario is sketched in Fig. 1b. A two-dimensional model of the melting region is justified because mantle flow is driven predominantly by plate-spreading and there is negligible plume-driven flow of mantle parallel to the axis (Section 2). We use the same thermal equation, lateral and upper boundary conditions as Bown and White (1995), as implemented by Walters et al. (2013). Model geometry and boundary conditions are illustrated in Fig. S3 (Supplementary Material). Plate-driven upwelling beneath a mid-ocean ridge is modelled as flow of an isoviscous fluid in a triangular corner with a wedge angle α = 70 • , within the range of values tested by Bown and White (1994), and a plate half spreading rate of U 1/2 = 10 km My −1 . This style of model has been extensively tested and applied to global compilations of oceanic crustal thickness in a range of tectonic settings White, 1994, 1995;Barton and White, 1997). Melting is calculated using the anhydrous parameterisation of Katz et al. (2003) with 15% modal clinopyroxene. The change of entropy on melting was set to 350 J kg −1 K −1 (Kojitani and Akaogi, 1997). Local instantaneous melt production rate, Γ , is calculated using Eq. (5) of Bown and White (1995). Both instantaneous melt extraction and melt emplacement focused at the spreading axis are assumed.
Temperature fluctuations at the basal boundary are imposed using where T 0 is the real mean temperature at the basal boundary, T is the amplitude of the temperature fluctuation, τ is the period of the temperature fluctuations and t is time. We specify τ = 5 Myr based on VSR observations (Jones et al., 2002b). T 0 is determined from the potential temperature at the basal boundary T P 0 , which is specified for each model. (T P 0 , T ) combinations were chosen to match the crustal thickness model in Fig. 3 and are listed in Table S2 (Supplementary Material). T P 0 varies across the study area from 1327 • C at the MAR/VSR-1 intersection to 1307 • C at the MAR/VST-2 intersection. This plume head temperature profile is slightly lower that estimated by  owing to the new crustal thickness model. T varies across the study area from 27 • C to 11 • C. The former figure agrees with previous estimates of 30 • and 25 • C for the average temperature change within the melting region at the northern end of the study area Poore et al., 2011). The corresponding full temperature variation at the base of the melting region varies from 54 • to 22 • C from north to south across the study area.
For adiabatic decompressional melting, the average composition C of the erupted melt is where c is the local instantaneous melt composition, X is the degree of melting and A is the area of the melting region. The assumption of adiabatic upwelling breaks down only in the top 5 km of the melting region, where the effect of conductive cooling is felt (Bown and White, 1994). About 80% of the total melting and much of the extraction of incompatible trace elements occurs in the adiabatic part of the melting region, so we consider Eq.
(2) a satisfactory approximation. To calculate c( X(z), T P 0 ), the mid-ocean ridge melting model was run to steady state temperature for a series of fixed (time-independent) T P 0 values between 1250 • and 1400 • C at 10 • C intervals. The horizontally averaged X(z) curve from each model was then used to calculate instantaneous melt compositions for various trace elements assuming fractional melting, using the forward modelling scheme of McKenzie and O'Nions (1991) and partition coefficients for aluminous lherzolite compiled by Gibson and Geist (2010). Our conclusions are not altered if trace element partition coefficients of McKenzie and O'Nions (1991) are used instead. We assume a homogeneous mantle source for simplicity in this first attempt at modelling VSRs. Compositional ranges for each element were gridded (Fig. S4, Sup-plementary Material). The c(z, T P 0 ) grids were then read back into the mid-ocean ridge melting model as a series of look-up tables. For time-dependent calculations, T P 0 was tracked across the calculation grid at each time-step, and the look-up tables were used to find c(x, z, t) for every point in the melting region and integrated to give C .

Results
A series of snapshots of evolving temperature within the melting region is shown in Fig. 8a. The steady-state temperature field has been subtracted from each snapshot to highlight temperature fluctuations. The sinusoidally varying basal temperature boundary condition generates thermal pulses that travel upward within the melting region, turn and spread laterally to become incorporated into mantle lithosphere. The boomerang shape of the thermal pulses arises because the upwelling rate is greatest in the centre of the melting region and decreases to zero at the base of the lithosphere. The depth to the base of the melting region moves up and down in response to the temperature fluctuations.
Predicted temporal evolution of crustal thickness, individual trace element concentrations and trace element ratios is plotted in Fig. 8. Melt production peaks correspond to times when a hot pulse lies about 10 km above the base of the melting region. Fluctuation in melt production is slightly skewed, with increases in melt production occurring more rapidly than decreases. This model prediction may explain some of the observed asymmetry in VSRs 1 and 3 visible in Fig. 8 of Poore et al. (2009). The younger flanks of the V-shaped ridges and older flanks of the troughs are observed to be steeper than the opposite flanks.
Compositions of individual elements vary at the 5 My periodicity of the imposed thermal fluctuation. Inverse correlation between composition and melt thickness is predicted. Associated maxima and minima vary in timing by up to 1 My. These phase differences reflect relationships between instantaneous melt composition, c, and depth (Fig. S4, Supplementary Material). For example, La and Yb both show strong (though opposite) variation in c near the base of the melting region, whereas c varies strongly at the top of the melting region for Sm. Therefore the Sm composition curve lags both the La and Yb curves. Phase difference between fluctuations in crustal thickness and trace element composition appears as closed curves on Fig. 8e. Stronger correlation between composition and crustal thickness is predicted for single elements than for element ratios.
Modelled compositional and crustal thickness relationships are compared with observations in Fig. 7. To generate an along-axis model prediction, composition-crustal thickness pairs were sampled at crustal thickness maxima of the VSR-1 and VSR-2 models and crustal thickness minima of the VST-1 and VST-2 models in Fig. 8. These four points were then interpolated to give model lines in Fig. 7. The melting model provides a good match to the observed relationship between crustal thickness and composition.
Amplitudes of observed and modelled compositional fluctuations are compared in Fig. 9. Amplitudes are expressed as minimum:maximum concentration. This presentation allows observations (recorded as concentrations) to be compared directly with model results (calculated as ratios of concentrations to source concentrations) without the need to specify the model homogeneous source composition. Details of the procedure used to make Fig. 9 and plots of model fit to observations are provided in Section C (Supplementary Material). The model provides a good match to across-axis compositional variation between VST-1 and VSR-2 for all representative single elements and element ratios. Alongaxis compositional variation between VST-1 and VSR-1 is generally well matched, except for ratios involving the most incompatible elements Nb/Y and La/Sm. Along-axis compositional variation between VST-1 and VSR-2 is less well matched, with discrepancies for ratios Nb/Y and Nb/La and single elements Zr, Sm, Y and Yb. The amplitude of the modelled compositional signal is compared with the amplitude of observed noise arising from fractional crystallization and incomplete magma mixing in Fig. 9b. Low signal:noise ratios for indicators Nb/Y, Nb/La and La/Sm probably explain why the model predictions lie slightly outside the uncertainty range of the along-axis data. Only the matches involving Zr, Sm, Y and Yb variation between VST-1 and VSR-2 may represent a real signal that is not accounted for in modelling, since signal:noise ratios for these elements are relatively high.
We conclude that the new melting model gives a generally good match to observed crustal thickness and trace element concentrations. Much of the VSR compositional signature can be explained by fluctuations in potential temperature of the mantle source.

VSR formation hypotheses
Although a thermal pulsing explanation of the VSR melting anomalies is only now becoming secure, the hypothesis of outward propagating melting anomalies had been generally accepted until Hey et al. (2010) and Benediktsdóttir et al. (2012) suggested that the VSRs result instead from rift propagation cycles. They showed that crustal accretion along the Reykjanes Ridge is asymmetric and postulated that the asymmetry arises from a series of rift relocations. They interpreted V-shaped ridge crests as pseudo-faults that mark past positions of southward-propagating rift tips. The geochemical data do not support this interpretation. Melt produced at propagating young rifts is known to be more enriched in trace elements than melt from equivalent established rifts because the average depth of melting beneath a propagating rift is deeper where the rift propagates into cooler, thicker lithosphere. The southern tip of the southward propagating Eastern Volcanic Zone on Iceland is a case in point (Walters et al., 2013). Similar enrichment in trace elements occurs near oceanic fracture zones (White et al., 1992). If the propagating rift scenario for the VSRs were correct then crustal thickness and trace element concentrations would be positively correlated, but observations show an inverse correlation. It therefore seems more likely that asymmetric crustal accretion along the Reykjanes Ridge occurs in addition to the thermal pulsing that generates VSR crustal thickness and geochemical signals. Asymmetric accretion appears to operate on a shorter timescale than VSR cycles.

Icelandic and Hawaiian plume fluxes
The new geochemical data show that VSR-forming melting anomalies propagate out to at least 1350 km from Iceland, well into the region where fracture zones obscure the topographic and gravitational expressions of the VSRs. This observation agrees with the plume head radius of 1800 km obtained from the thermally supported plume swell (Fig. 2a). The map-view geometry of the VSRs can be used to infer propagation rates of the thermal pulses within the asthenosphere, and hence plume flux (Vogt, 1971;White and Lovell, 1997). Propagation rates measured from VSR topography are 87 to 255 km yr −1 , which correspond to plume volume fluxes of 28 to 90 km 3 yr −1 if the plume head layer is 100 km thick (Poore et al., 2009). The four VSR geometry models shown on Fig. 4 (bottom right panel) have fluxes of 50 ± 13, 49 ± 14, 49 ± 7 and 47 ± 5 km 3 yr −1 (Section D, Supplementary Information). We adopt the conservative flux value of 49 ± 14 km 3 yr −1 which covers all four models, as well as the mean value of Poore et al. (2009). The gravity field has lower spatial resolution than the topography data used by Poore et al. (2009). Nevertheless, gravity-derived VSR flux estimates are perhaps more reliable than those derived from topography because the lower spatial resolution helps to filter out the asymmetric crustal accretion effects noted by Hey et al. (2010) and Benediktsdóttir et al. (2012).
Another class of plume flux estimates comes from models of the entire plume head. Estimates in this class have consistently been smaller than VSR-based estimates. Sleep (1990) estimated a volume flux of 0.3 km 3 yr −1 using a flux balance method. Ribe and Delattre (1998) estimated ∼1 km 3 yr −1 using a lubrication theory model. Ito (2001) estimated 5.7 km 3 yr −1 using a 3D numerical model. These plume head models assume along-axis interaction distances of 400, 460 and 500 km of respectively, all considerably smaller than the radius of 1800 km estimated from the plume swell and supported by joint geochemical and geophysical VSR observations. Some of the small plume radius estimates were influenced by the region of strong trace element and isotopic enrichment surrounding Iceland (Fig. 2b). However, the edge of the strongly enriched region does not correspond to the full radius of the plume head (Taylor et al., 1997;Murton et al., 2002;Shorttle et al., 2010). Since plume flux is proportional to the square of the plume head radius, underestimates of plume-head radius are expected to have a significant effect on plume flux estimates. Sim- ple scaling of the volume flux estimates listed above to a radius of 1800 km gives revised estimates of 6, 15 and 74 km 3 yr −1 . Hence, incorrect identification of the edge of the plume head goes a long way to explaining the order-of-magnitude difference between VSRderived and plume head model-derived flux estimates. Until whole plume-head models of the correct radial size are developed, VSRderived plume flux estimates should be preferred.
U-series isotopic disequilibria in young lavas provide a measure of mantle upwelling velocity within the deep part of the melting region. Koornneef et al. (2012) estimate an upwelling rate of 112 mm yr −1 within the Iceland Plume conduit. Since the plume conduit has a radius of about 150 km (Jones and Maclennan, 2005;Shorttle et al., 2010), an estimate of volume flux through the melting region is 7.9 km 3 yr −1 . Jones and Maclennan (2005) provide another estimate of the melt flux from deep, plume-driven upwelling of ∼0.1 km 3 yr −1 based on seismic crustal thickness measurements at Iceland. Since ∼2% melting occurs in this deep part of the melting region (Maclennan et al., 2001), the flux of mantle through deep melting region is ∼5 km 3 yr −1 , which is similar to the value of 7.9 km 3 yr −1 obtained independently from U-series isotopes. Volume flux through the melting region is expected to be less than the total plume flux because not all mantle welling up the plume conduit melts. The base of the melting region is shallower than the base of the asthenosphere and the diameter of the melting region is narrower than the diameter of the plume conduit (Fig. 15 of Watson and McKenzie, 1991). Considering ranges of 49 ± 14 km 3 yr −1 for the total flux and 5.0-7.9 km 3 yr −1 for the flux that melts, the proportion of plume mantle processed by melting at the top of the conduit is 8-22%. Most of this processing is by small degrees of deep melting, in the presence of garnet. Such prior processing of the mantle that arrives beneath our study area can probably explain the depletion in light and medium REE noted by Shorttle et al. (2010).  Kokfelt et al. (2003); circles, Stracke et al. (2003Stracke et al. ( , 2006; squares, Koornneef et al. (2012). Hawaiian data (open symbols): circles, Cohen and O'Nions (1993); squares, Sims et al. (1995Sims et al. ( , 1999; diamonds, Pickett and Murrell (1997). All activity ratios adjusted to decay constants in Koornneef et al. (2012). Lines are ±1 standard deviation from a linear regression: solid lines, Iceland; dashed lines, Hawaii. Distances calculated from the plume centres determined from both geophysical and geochemical data (Watson and McKenzie, 1991;Shorttle et al., 2010).
The Hawaiian Plume has been considered the most vigorous mantle plume (Sleep, 1990;Bourdon et al., 2006). Models of the Hawaiian plume swell and melt production rate give volume flux estimates of 9.5, 7.9 and 3.7 km 3 yr −1 (Sleep, 1990;Watson and McKenzie, 1991;Ribe and Christensen, 1994). These fluxes are significantly less than our estimate of 49 ± 14 km 3 yr −1 , and suggest that the Iceland Plume is in fact the most vigorous plume in Earth's mantle at present. This conclusion is contentious because U-series disequilibrium measurements at Iceland have been quoted as being greater than those at Hawaii, which would suggest that the average upwelling rate at the base of the melting region is greater beneath Hawaii . We have assembled a dataset of U-series measurements from Iceland and Hawaii to investigate this issue (Fig. 10). All data series show a component of increase in activity ratio with distance from the plume centre, which we quantify using linear regressions. We interpret these activity ratio-distance relationships as the signature of plume-driven upwelling, following Kokfelt et al. (2003), Bourdon et al. (2006) and Koornneef et al. (2012). There is no significant difference between activity ratio-distance relationships at Iceland and Hawaii for ( 231 Pa/ 235 U) data (F-test, 99.5% probability) or ( 230 Th/ 238 U) data (98.0% probability). Higher average values of these activity ratios quoted by Bourdon et al. (2006) for Iceland in comparison with Hawaii reflect the larger number of Icelandic samples at radial distances of >200 km, outside the influence of significant plume-driven upwelling. Dynamic melting calculations show that ( 231 Pa/ 235 U) and ( 230 Th/ 238 U) are relatively insensitive to variation in mantle upwelling rates above 100 mm yr −1 , equivalent to volume fluxes above ∼10 km 3 yr −1 Stracke et al., 2006). Available U-series data are therefore consistent with our conclusion that the Iceland Plume has a greater flux than the Hawaiian Plume. We consider that 49 ± 14 km 3 yr −1 is the best current estimate for the volume flux of the Iceland Plume because it has the best potential to reconcile flux estimates from VSR geometry, plume head models and U-series geochemistry.

Conclusions
Along the Reykjanes (Mid Atlantic) Ridge south of Iceland, spatial variations in incompatible trace element concentrations and ratios correlate with V-Shaped Ridge (VSR) records from bathymetry and gravity. Two complete cycles of variation are observed for the first time. Chemical variation between a VSR crest and the adjacent trough on the same plate spreading flowline is also observed for the first time. Compositional variations associated with the VSRs cannot be explained by fractional crystallization because there is no corresponding variation in Mg number. Enrichment in incompatible trace elements is inversely correlated with estimates of VSR crustal thickness from wide-angle seismic and gravity data. This observation implies a positive correlation between average melt fraction and crustal thickness, so that temperature variation in the source mantle is important in controlling crustal thickness variations.
A model of mid-ocean ridge melting along a two-dimensional transect parallel to the plate spreading direction has been developed to investigate VSR melt thickness and composition. Sinusoidal temperature variation is imposed at the base of the melting region, and the trace element composition of the melt is calculated by accounting for time-dependent temperature variations throughout the melting region. The new model gives a good match to the observed inverse correlation between crustal thickness and trace element concentrations. Observed amplitudes of VSR-related compositional variation are mostly well matched by the model.
The new observations and models suggest that VSRs are generated when hotter and cooler blobs are carried up the plume stem beneath SE Iceland and spread radially outward within the asthenosphere, before being drawn upward into the melting region when passing beneath the Mid Atlantic Ridge. The observed inverse correlation between crustal thickness and trace element concentrations is not compatible with a recent hypothesis that VSRs are formed by propagating mid-ocean ridge segment offsets. Peak-to-peak temperature variation associated with plume pulsing varies from 54 • C at a distance of 730 km from the plume centre to 22 • at a distance of 1280 km. These temperature estimates can now be used to produce detailed models of vertical motion of the Greenland-Scotland Ridge, which will in turn allow more detailed analysis of the link between plume pulsing and Northern Component Water flux into the global ocean through the Neogene.
The VSR compositional signal is observed out to at least 1350 km from Iceland, in agreement with a plume head radius of 1800 km estimated from the plume swell. Plume volume flux estimated from the map-view geometry of the VSRs is 49 ± 14 km 3 yr −1 . Some previous flux estimates an order of magnitude lower are not reliable because they assumed an unrealistically small plume head. 8-22% of plume mantle is processed by melting at the top of the conduit. The Iceland Plume is the most vigorous plume in Earth's mantle at present. U-series isotope data from Iceland and Hawaii are consistent with this hypothesis.