Rate of Melt Ascent beneath Iceland from the Magmatic Response to Deglaciation

TL;DR

This study models melt ascent velocity beneath Iceland during deglaciation, revealing that a finite melt transport speed explains observed REE variations and volcanic eruption timing better than instantaneous models.

Contribution

It introduces a numerical model incorporating finite melt ascent velocity, improving understanding of volcanic responses and REE variations during deglaciation.

Findings

Melt ascent velocity during deglaciation is estimated to be around 100 m/yr.

Finite melt transport explains large REE variations in Icelandic lavas.

Previous models with instantaneous transport underestimate REE depletion.

Abstract

Observations of the time lag between the last deglaciation and a surge in volcanic activity in Iceland constrain the average melt ascent velocity to be $\geq50$ $\mathrm{m/yr}$. Although existing theoretical work has explained why the surge in eruption rates increased $5$-$30$ fold from the steady-state rates during the last deglaciation, they cannot account for large variations of Rare Earth Element (REE) concentrations in the Icelandic lavas. Lavas erupted during the last deglaciation are depleted in REEs by up to $70\%$; whereas, existing models, which assume instantaneous melt transport, can only produce at most $20\%$ depletion. Here, we develop a numerical model with finite melt ascent velocity and show that the variations of REEs are strongly dependent on the melt ascent velocity. When the average melt ascent velocity is $100$ $\mathrm{m/yr}$, the variation of $\mathrm{La}$ calculated by our model is comparable to that of the observations. In contrast, when the melt ascent velocity is $1,000$ $\mathrm{m/yr}$ or above, the model variation of $\mathrm{La}$ becomes significantly lower than observed, which explains why previous models with instantaneous melt transport did not reproduce the large variations. We provide the first model that takes account of the diachronous response of volcanism to deglaciation. We show by comparing our model calculations of the relative volumes of different eruption types (subglacial, finiglacial and postglacial) and the timing of the bursts in volcanic eruptions with the observations across different volcanic zones that the Icelandic average melt ascent velocity during the last deglaciation is likely to be $\sim100$ $\mathrm{m/yr}$.

Tables

Table 1. Table 1 : Parameter values for calculations.

Parameter	Meaning	Value	Dimensions
$D^{\mathrm{La}}$	$\mathrm{La}$ partition coefficient	$0.010$	$1$
$-{\left(\partial F/\partial P\right)}_S$	isentropic melt productivity	10	$\mathrm{wt}\%/\mathrm{GPa}$
$g$	gravitational acceleration	$9.82$	$m/s^2$
$P_{sol}$	solidus pressure	$3.5$	$\mathrm{GPa}$
$U_0$	plate half-spreading rate	$10$	$\mathrm{mm}/\mathrm{yr}$
$z_c$	crustal thickness	$20$	$\mathrm{km}$
$\alpha$	ridge angle	$45$	$\mathrm{deg}$
$\eta$	mantle viscosity	$8.0\times {10}^{18}$	$\mathrm{Pa}s$
${\rho }_i$	density of ice	$900$	$\mathrm{kg}/m^3$
${\rho }_l$	density of melt	$2900$	$\mathrm{kg}/m^3$
${\rho }_s$	density of solid mantle	$3300$	$\mathrm{kg}/m^3$
${\tau }_B$	yield stress of ice	$100$	$\mathrm{kPa}$

Table 2. Table 2 : Timeline of the model ice sheet.

Geological Period	Time (kyrBP)	Ice Radius (km)	Average Ice Thickness (km)
Last Glacial Maximum	before $23.0$	$300.0$ constant	$1.39$ constant
Offshore Deglaciation	$23.0\to 17.0$	$300.0\to 180.0$	$1.39\to 1.08$
Pre Bølling-Allerød	$17.0\to 15.0$	$180.0$ constant	$1.08$ constant
Bølling-Allerød	$15.0\to 13.8$	$180.0\to 160.0$	$1.08\to 1.02$
Younger Dryas	$13.8\to 11.7$	$160.0$ constant	$1.02$ constant
Early Holocene	$11.7\to 10.5$	$160.0\to 45.0$	$1.02\to 0.54$
Holocene	$10.5\to$ now	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}45.0}$ constant	$0.54$ constant

Table 3. Table 3 : Model distances of the zones from ice center.

Zone	Ranges (km)
WVZN	$0$–$70$
WVZS	$70$–$180$
NNVZ	$120$–$180$
REYK	$180$–$250$

Equations

$$h(r,t)=\{\begin{array}{ll}{h}_m(t)\sqrt{1-\frac{r}{r_m(t)}},& 0\le r\le r_m(t),\\ 0,& \text{otherwise,}\end{array}$$

$$h_m(t)$$

$$r_m(t)$$

$$\frac{DF}{Dt}={\left(\frac{\partial F}{\partial P}\right)}_S\frac{DP}{Dt}$$

$$\Gamma (\mathbf{x},t)$$

$$={\rho }_s{\left(\frac{\partial F}{\partial P}\right)}_S\frac{DP}{Dt}.$$

$$F(P)=-{\left(\frac{\partial F}{\partial P}\right)}_S\left(P_{sol}-P\right).$$

$$\dot{M}(t)$$

$$\dot{V}(t)$$

$$=\frac{1}{{\rho }_l}{\int }_{\mathcal{V}}\Gamma \left(\mathbf{x},t-\frac{|z|}{v_t}\right)dV.$$

$$\frac{c_l^i}{c_{s0}^i}=\frac{1}{D^i}{\left(1-F\right)}^{\frac{1}{D^i}-1}$$

$$\dot{M_i}(t)$$

$$\dot{V_i}(t)$$

$$=\frac{1}{{\rho }_l}{\int }_{\mathcal{V}}{c}_l^i\Gamma \left(\mathbf{x},t-\frac{|z|}{v_t}\right)dV.$$

$${\bar{c}}_l^i(t)=\frac{{\dot{M}}_i(t)}{\dot{M}(t)}=\frac{{\dot{V}}_i(t)}{\dot{V}(t)}.$$

$$v_x(x,y,z)$$

$$v_y(x,y,z)$$

$$\text{and}{v}_z(x,y,z)$$

$$B$$

$$v_r(r,z,t)=$$

$$v_{\theta }(r,z,t)=$$

$$v_z(r,z,t)=$$

$$P(r,z,t)=$$

$$\text{and}\frac{DP}{Dt}(r,z,t)=$$

$$\tau \equiv$$

$$\tilde{w}(k,t)=$$

$${\mathcal{H}}_n[f](k)=$$

$${\mathcal{H}}_n^{-1}[F](r)=$$

$${\mathcal{H}}_0[\dot{h}](k,t)=\frac{3\sqrt{2}}{32}\dot{V}(t)\left[{\mathrm{J}}_{-\frac{1}{4}}\left(\frac{k{r}_m(t)}{2}\right){\mathrm{J}}_{\frac{1}{4}}\left(\frac{k{r}_m(t)}{2}\right)+{\mathrm{J}}_{-\frac{3}{4}}\left(\frac{k{r}_m(t)}{2}\right){\mathrm{J}}_{\frac{3}{4}}\left(\frac{k{r}_m(t)}{2}\right)\right].$$

Full text

Key Points:

• •

We model the magmatic response to the last deglaciation in Iceland with finite melt ascent velocity.

• •

The model results compared with observations suggest that the melt ascent velocity is likely to be around $100$ $\mathrm{m/yr}$.

1 Introduction

Iceland is located where a mantle plume meets the Mid-Atlantic Ridge (White \BOthers., \APACyear1992). The mantle plume and the spreading center are responsible for the upwelling of the mantle underneath Iceland, which induces decompression melting in the upper mantle (McKenzie \BBA Bickle, \APACyear1988). This decompression melting produces magma that supplies the production of the Icelandic crust and volcanic eruptions. Geological observations indicate that eruption rates in the volcanic zones of Iceland were significantly elevated during a burst of activity that took place after the end of the last major deglaciation (Sigvaldason \BOthers., \APACyear1992; Slater \BOthers., \APACyear1998; Maclennan \BOthers., \APACyear2002; Sinton \BOthers., \APACyear2005; Maclennan, \APACyear2008; Eason \BOthers., \APACyear2015). This period of high productivity, perhaps $30$ times or more, may have started about $15$ $\mathrm{kyrBP}$ and ended before $9$ $\mathrm{kyrBP}$.

The cause of the surge in eruption rates has been examined by Jull \BBA McKenzie (\APACyear1996). In their model, post-glacial rebound induced by the last deglaciation increases the rate of pressure drop in the upper mantle by up to $50$-fold from the steady-state value. This increase in the decompression rate significantly increases the melting rate in the upper mantle, which leads to greater melt supply.

Jull \BBA McKenzie (\APACyear1996) also showed that the decompression rate due to post-glacial rebound has its maximum value at the surface and decays exponentially with depth. This means that the additional melt production in the mantle during the deglaciation occurs mostly at shallow depths. The melts generated at these depths are depleted in Rare Earth Elements (REEs). These additional melts produced during the deglaciation will therefore dilute the concentrations of REEs in the aggregated melts. By assuming that the melt transport is instantaneous, Jull \BBA McKenzie (\APACyear1996) calculated that the REE concentrations in the melts decrease by around $20\%$ during the last deglaciation compared to melts generated at other times when the ice-load is thought to have been close to steady-state. However, geological observations indicate that lavas erupted during the surge in volcanic eruption rates are depleted in REEs by up to $\approx 70\%$ (Slater \BOthers., \APACyear1998; Maclennan \BOthers., \APACyear2002; Sinton \BOthers., \APACyear2005; Maclennan, \APACyear2008; Eason \BOthers., \APACyear2015), which is significantly higher than that calculated by Jull \BBA McKenzie (\APACyear1996).

Slater \BOthers. (\APACyear1998) attempted to account for this mismatch by developing an inverse model similar to that of McKenzie \BBA O’Nions (\APACyear1991), which used the observed variations of REE concentrations to constrain the melt productivity function. Slater \BOthers. (\APACyear1998) showed that there exists a melt productivity function that matches the Jull \BBA McKenzie (\APACyear1996) theoretical variations of the REE concentrations with the geological data. However, in order for such a melt productivity function to exist, model parameters had to be modified including reducing the initial ice sheet radius to $90$ $\mathrm{km}$, which is significantly smaller than the likely radius of the ice sheet as inferred from observational studies (Sigmundsson, \APACyear1991; Hubbard \BOthers., \APACyear2006; Licciardi \BOthers., \APACyear2007; Pétursson \BOthers., \APACyear2015; Patton \BOthers., \APACyear2017).

Maclennan \BOthers. (\APACyear2002) previously used the relative timing of the last deglaciation and the surge in volcanic eruption rates in the Northern Volcanic Zone (NVZ) of Iceland to estimate that the melt ascent velocity is at least $50$ $\mathrm{m/yr}$. The Icelandic melt ascent velocity during the mid-Holocene has also been estimated from a time lag of $\approx 600$ $\mathrm{yr}$ (Swindles \BOthers., \APACyear2017) obtained from the cross-correlation between the Icelandic volcanic eruption rates and the change in the atmospheric circulation pattern indicated by sodium concentrations in Greenland Ice Sheet Project 2 (GISP2). This $\approx 600$ $\mathrm{yr}$ time lag gives an estimated Icelandic melt ascent velocity of $\approx 50$–$100$ $\mathrm{m/yr}$ during the mid-Holocene.

Here, we investigate how the melt ascent velocity in the mantle and the crust influences the variations of REE concentrations. By incorporating a finite rate melt transport model into the model of Jull \BBA McKenzie (\APACyear1996), we show that variations of REE concentrations depend significantly on the melt ascent velocity. With an appropriate melt ascent velocity, our model demonstrates that the model variations of REE concentrations can be matched with those observed geologically. While Jull \BBA McKenzie (\APACyear1996) used a very simple ice-load model with a constant ice radius, our model combines an ice-load history with melt generation and transport and therefore enables prediction of the volumes and compositions of melt that are erupted either subglacially or in ice-free settings. This feature of the model allows for direct comparison with geological observations, which use edifice geomorphology and volcanic facies analysis to determine whether an eruption is subglacial or not. Therefore, not only does this work help us understand how melt transport affects REE concentrations and eruption types in different places on Iceland, but also it can be useful as a tool to constrain the melt transport rate.

In the next section, we will begin by considering how the mantle flow responds to deglaciation and use this response to calculate the melting rate and the compositions of the melts generated. We then transport the melts produced at depth to the surface to calculate the eruption rate together with the composition of the erupted lavas. Finally, we compare our numerical results to the observational data in order to constrain the melt ascent velocity.

2 Model

We follow the modelling of Jull \BBA McKenzie (\APACyear1996) with the following key differences:

While Jull \BBA McKenzie (\APACyear1996) used an ice-load with a constant radius, our ice sheet behaves like a gravity current with time-dependent radius and thickness. 2. 2.

Our ice-load input consists of multiple deglaciation stages beginning at $23.0$ to $10.5$ $\mathrm{kyrBP}$ designed to capture key features of ice-sheet reconstructions. 3. 3.

We neglect the elastic response of the solid mantle. 4. 4.

We assume finite melt ascent velocity.

Numerical parameters used as inputs into the model are listed in Table 1.

Numerical values of the plate half-spreading rate ($U_{0}$), crustal thickness ($z_{c}$), ridge angle ($\alpha$) and mantle viscosity ($\eta$) are the same as in Jull \BBA McKenzie (\APACyear1996). $z_{c}=20$ $\mathrm{km}$ is at the lower bound of the Darbyshire \BOthers. (\APACyear2000) estimates ($20$–$37$ $\mathrm{km}$) because our study areas are relatively far ($>100$ $\mathrm{km}$) from the mantle plume center. Numerical values of $U_{0}=10$ $\mathrm{mm/yr}$ and $\eta=8.0\times 10^{18}$ $\mathrm{Pa\,\,s}$ are similar to the Árnadóttir \BOthers. (\APACyear2009) estimates. The density of ice ($\rho_{i}$) is in the range of $830$–$917$ $\mathrm{kg/m^{3}}$ in Paterson (\APACyear1994). The densities of melt ($\rho_{l}$) and of solid mantle ($\rho_{s}$) follow Katz \BOthers. (\APACyear2003). Sources of the remaining numerical parameters will be mentioned later.

2.1 Glacial Load

Due to limited geological records of the ice sheet, it is not straight-forward to reconstruct the details of the shape of the ice-sheet during the last deglaciation (Hubbard \BOthers., \APACyear2006; Patton \BOthers., \APACyear2017). In Jull \BBA McKenzie (\APACyear1996), the ice sheet was assumed to have axisymmetric parabolic shape with a constant radius of $180$ $\mathrm{km}$. Here, we modify the ice sheet to be an axisymmetric viscous gravity current (Paterson, \APACyear1994) with glacier terminus retreating during the deglaciation. This is a more reasonable representation of the actual ice sheet, allowing the spatial variations of the volcanic response to be examined more accurately.

In an axisymmetric gravity current ice model (Huppert, \APACyear1982; Paterson, \APACyear1994), the differential thickness of ice along the radial direction induces a radially-inward basal shear stress. At the yield strength limit of ice, the basal shear stress is uniform and is equal to the yield stress of ice. When the stress exceeds the yield strength limit, ice deformation and sliding will occur. These processes will re-adjust the ice sheet shape until it returns back to within its yield strength limit. We assume here that the time scale of the ice deformation and sliding (when the stress exceeds the yield strength) is short compared to the time scale of the deglaciation. That is, the glacier is assumed to always be within its yield strength limit and the thickness of ice $h(r,t)$ as a function of radial distance $r$ and time $t$ follows (Huppert, \APACyear1982; Paterson, \APACyear1994)

[TABLE]

where

[TABLE]

$h_{m}(t)$ is the thickness of ice at the center, $r_{m}(t)$ is the radial extend of ice, $V(t)$ is the volume of ice, $\rho_{i}$ is the density of ice and $\tau_{B}$ is the yield stress of ice.

The numerical value of $\tau_{B}=100$ $\mathrm{kPa}$ we use here (Paterson, \APACyear1994) gives ice sheet dimensions that closely resemble that of the Late Weichselian Icelandic ice sheet (Sigmundsson, \APACyear1991; Hubbard \BOthers., \APACyear2006; Licciardi \BOthers., \APACyear2007; Pétursson \BOthers., \APACyear2015; Patton \BOthers., \APACyear2017) and can reproduce the ice radius of $180$ $\mathrm{km}$ together with $2$ $\mathrm{km}$ ice thickness used previously in Jull \BBA McKenzie (\APACyear1996).

The time evolution of the ice coverage as an input into our model follows approximately that of Patton \BOthers. (\APACyear2017). We set the input deglaciation to consist of three stages during which the ice volume decreases linearly with time and the ice volume stays constant during two intermissions between these three deglaciation stages. Time $t=0$ in the model corresponds to the present (AD 1950).

We set the initial ice load to have a radius of $300$ $\mathrm{km}$ covering the whole of Iceland and most of the continental shelf. The ice volume is held constant until $t=-23.0$ $\mathrm{kyr}$ when the first stage of deglaciation (we refer to as the Offshore Deglaciation) begins. The Offshore Deglaciation terminates at the shoreline with ice radius of $180$ $\mathrm{km}$ at time $t=-17.0$ $\mathrm{kyr}$ followed by a pause of $2.0$ $\mathrm{kyr}$. Next, the second stage of deglaciation (Bølling-Allerød) proceeds from time $t=$ $-15.0$ to $-13.8$ $\mathrm{kyr}$ during which the ice radius decreases from $180$ $\mathrm{km}$ to $160$ $\mathrm{km}$. Then, the deglaciation pauses for $2.1$ $\mathrm{kyr}$, corresponding approximately to the Younger-Dryas. The final stage of deglaciation (Early Holocene) takes place between time $t=$ $-11.7$ and $-10.5$ $\mathrm{kyr}$ with the ice sheet retreating from radius of $160$ $\mathrm{km}$ to $45$ $\mathrm{km}$, which is approximately the current size of the Vatnajökull ice sheet. The timeline of the modelling-input deglaciation is comparable to the last deglaciation in Iceland (Sigmundsson, \APACyear1991; Maclennan \BOthers., \APACyear2002; Hubbard \BOthers., \APACyear2006; Licciardi \BOthers., \APACyear2007; Pétursson \BOthers., \APACyear2015; Patton \BOthers., \APACyear2017) and is summarised in Figure 1, Table 2 and the Supporting Information.

2.2 Mantle Flow

Similar to the assumption made by Jull \BBA McKenzie (\APACyear1996), in steady state, the spreading ridge induces passive upwelling of the mantle, which we assume to follow corner flow (Batchelor, \APACyear2000; Spiegelman \BBA McKenzie, \APACyear1987). Active upwelling induced by the mantle plume can also increase the melt production rate. However, the geological data in our study come from regions that are at least $\sim 100$ $\mathrm{km}$ away from the plume center. We therefore assume that the active upwelling is insignificant here.

The glacial load on the surface affects the pressure in the mantle underneath. During deglaciation, the surface load drops, which leads to an increased mantle decompression melting rate from that induced by the steady state passive corner flow (Figure 2).

To calculate the effect of deglaciation on mantle melting rate, we first note that the Maxwell relaxation time ($\tau_{M}=\eta/\mu\approx 10^{1}$ $\mathrm{yrs}$) is much shorter than the viscous characteristic time ($\tau_{v}=2\eta k/\rho_{s}g\approx 10^{3}$ $\mathrm{yrs}$) where $\mu$ is the elastic modulus, $k$ is wave number and other variables are as defined in Table 1. This means that the elastic deformation in the viscoelastic mantle model used in Jull \BBA McKenzie (\APACyear1996) is negligible and the deformation in the mantle is dominated by the viscous response. We therefore model the mantle as a viscous half-space incompressible fluid and the elastic thickness of the Icelandic lithosphere is assumed to be negligible. When using the same modelling inputs as in Jull \BBA McKenzie (\APACyear1996), our numerical model yields the same results as those in Jull \BBA McKenzie (\APACyear1996), which verifies our assumption that the elastic deformation is insignificant.

The boundary conditions on the surface of the half-space mantle are that the normal stress is equal to the pressure from the weight of the ice load and that the shear stress is negligible compared to the normal stress. We obtain semi-analytical solutions to the mantle flow in cylindrical coordinates in response to the glacial load as provided in Appendix A together with the corner flow solutions.

2.3 Decompression Melting in the Mantle

Mantle upwelling is sufficiently fast that the heat loss due to conduction is negligible. Therefore, the decompression melting in the mantle is adiabatic and the mantle melting rate $DF/Dt$ can be calculated by (Jull \BBA McKenzie, \APACyear1996)

[TABLE]

where $F$ is the degree of melting by mass fraction relative to the initial mass of the solid mantle, $\left(\partial F/\partial P\right)_{S}$ is the isentropic melt productivity of the mantle and $\frac{D}{Dt}$ is the convective derivative following the solid mantle trajectories.

The rate of mass production of melt per unit volume as a function of space and time is assumed to follow

[TABLE]

The isentropic melt productivity depends on several factors including the composition of the mantle, temperature and pressure (McKenzie, \APACyear1984). In numerical calculations, using different melt productivity functions will result in different profiles of depth-dependent mantle melting rate, and different eruptive REE concentrations (Slater \BOthers., \APACyear1998). To investigate the effect of magma transport solely without the effect of melt productivity function on the eruptive REE concentrations, we use a constant isentropic melt productivity (Table 1) and the degree of melting as a function of pressure follows a simple linear relation

[TABLE]

Our choice of solidus pressure and melt productivity (in Table 1) gives a melt productivity function that closely resembles that obtained from the melt parametrisation of Katz \BOthers. (\APACyear2003) at $1500$ $\mathrm{{}^{o}C}$ potential temperature.

Sims \BOthers. (\APACyear2013) have shown that the temporal variability of isotope ratios in lavas erupted during the last deglaciation in northern Iceland provide evidence for a lithologically heterogeneous mantle source beneath Iceland. We investigate the effect of mantle heterogeneities by comparing our simple homogeneous mantle model results to that of the pMELTS modelling (Ghiorso \BOthers., \APACyear2002; Smith \BBA Asimow, \APACyear2005) of a bi-lithological mantle as used in Rudge \BOthers. (\APACyear2013). We show these results in Supporting Information that both mantle models yield the same conclusions for the rate of melt ascent. Our model is not very sensitive to the mantle heterogeneities because the model calculations do not involve isotopic composition.

Melts generated in the mantle have to be transported to the surface before they erupt. We assume that the effects of finite melt transport rate can be approximated by sampling the melt production rate field (equation (3)) with a time-lagged sampler. To the leading order, we assume that the vertical component of the melt velocity is constant $=v_{t}$. In this case, the time taken for melt produced at location $(x,y,z)$ in the mantle to ascend to the surface is $\Delta t=|z|/v_{t}$, where $|z|=-z$ ($\because z<0$ below the Earth’s surface). That is, melt that reaches the surface at time $t$ is assumed to have been produced at time $t^{\prime}=t-|z|/v_{t}$ in the past. Therefore, the total mass flux of melt supply to the crust at time $t$ is

[TABLE]

which is the integral of all the instantaneous melts produced in the melting region $\mathcal{V}$; however, the melts added from depth $|z|$ are assumed to have been produced at time $t^{\prime}=t-|z|/v_{t}$ in the past.

The total volume flux $\dot{V}$ of melt supply to the crust at time $t$ can be calculated from the mass flux:

[TABLE]

2.4 REE Concentrations

We simplify the model by assuming that the concentration $c^{i}_{l}$ of a highly incompatible element $i$ with partition coefficient $D^{i}$ in the instantaneous melt can be calculated based on modal fractional melting (Shaw, \APACyear1970)

[TABLE]

where $c^{i}_{s0}$ is the concentration of the element in the initial source.

Equation (7) gives the instantaneous concentration as a function of the degree of melting $c^{i}_{l}=c^{i}_{l}(F)$. The degree of melting as a function of pressure $F=F(P)$ is known from equation (4) and the pressure as a function of position and time $P=P(\boldsymbol{x},t)$ is known from equation (2). We can therefore combine these equations to calculate at any location in the mantle at any time the instantaneous concentration $c^{i}_{l}=c^{i}_{l}(\boldsymbol{x},t)$ in the melt generated. The bulk partition coefficient of $\mathrm{La}$ (Table 1) is assumed to follow that in Workman \BBA Hart (\APACyear2005).

This very simplified melting modelling of $\mathrm{La}$ gives results that are not significantly different from those (shown in Supporting Information) obtained from a more elaborate model of mantle melting used in Rudge \BOthers. (\APACyear2013) because highly incompatible elements (such as $\mathrm{La}$) partition into melts almost completely near the solidus intersection in the garnet field.

Given the concentration (by mass) $c^{i}_{l}$ of a trace element $i$ in the instantaneous melt as a function of space and time, the total mass flux of the trace element $i$ in the melt supply to the crust is

[TABLE]

where $c^{i}_{l}$ is calculated at point $(\boldsymbol{x},t-\frac{|z|}{v_{t}})$.

Similar to the volume flux of the whole melt defined in equation (6), we define the total “volume” flux of a trace element $i$ in the melt supply to the crust as

[TABLE]

Following these definitions, the mean concentration of the element $i$ in the melt supply to the crust at time $t$ is

[TABLE]

3 Results and Discussion

3.1 Decompression Melting and Eruption Rates

The numerical methods we use for the calculations are discussed in Appendix B. Figure 3 illustrates snapshots of the decompression rate in the mantle from the model when the ice load history follows the timeline given in Section 2.1.

While the Jull \BBA McKenzie (\APACyear1996) model with constant radius of ice-load predicted that the region of maximum decompression rate is always below the center of the ice sheet (their Figure 3), our model with variable ice radius predicts that this region is below the glacier terminus and is moving radially as the ice retreats. The glacially induced decompression causes the spatially dependent mantle melting rate underneath Iceland to increase from its steady state value by several fold during the deglaciation. These extra melts then transport to the surface, causing an increase in volcanic eruption rates.

The time delay between the surge of mantle melting and the surge of volcanic eruptions depends on the melt transport speed and also on how long melts reside in crustal chambers before they erupt. Figure 4a shows the melt supply rates to crustal chambers predicted by our model from different input values of melt ascent velocity by integrating equation (6) in the melting region underneath Iceland along the ridge axis from $45$ to $270$ $\mathrm{km}$ from the ice center, taking into account the time delay due to finite melt ascent velocity. The graph demonstrates that if melt transport were almost instantaneous, the surge in the melt supply rate (red curve) would respond almost immediately after the deglaciation period (grey shaded area). Whereas, with slower melt transport, the surge in the melt supply rate will be delayed from the deglaciation period. At lower rates of melt transport, the shape of the melt supply rate curve will be more stretched in time because melts produced at the same time at different depths will arrive at crustal chambers at different times.

Note that the area under the curve over the whole time interval shown in the graph is independent of the melt ascent velocity. This is because, by the conservation of mass, the total melt supply is equal to the total melt produced regardless of how fast the melt is transported.

Before melts erupt on the surface, their compositions can be modified in the crustal chambers. We assume that the amount of melts accommodated in a chamber is constant. By conservation of mass, this implies that the total mass flux into is equal to the total mass flux out of the chamber. Therefore, the eruption rate is equal to the rate of melts entering the chamber (Figure 4a). However, the mass flux of each individual component do not need to follow this rule. Mixing and crystallization processes can modify the concentrations of REEs. We will discuss these two processes together with the remaining plots in Figure 4 later in Section 3.5.

3.2 Eruptive Locations

The ages and volumes of eruptions from the last glacial and present postglacial are compiled using published maps and age estimates. The principal sources of information for the Northern Volcanic Zone (NVZ) are Sæmundsson (\APACyear1991) and Sæmundsson \BOthers. (\APACyear2012). For the Western Volcanic Zone (WVZ) and Reykjanes Peninsula (REYK), the maps of Sinton \BOthers. (\APACyear2005), Eason \BOthers. (\APACyear2015) and Sæmundsson \BOthers. (\APACyear2016) are used. Acknowledgements section lists all the sources of rock sample dataset we use in this work.

Geological mapping, geomorphology and interpretation of the volcanic lithologies have been used to determine the eruptive facies: whether it is subglacial, finiglacial or postglacial. Finiglacial means that there is evidence of thin or recently disappeared ice when the eruption unit was being formed. Finiglacial units are likely to have formed when the glacier terminus was sweeping through the eruptive area during the glacial retreat.

Tephrochronology provides bounds on eruption ages for the postglacial events, meaning that the age constraints are expressed with a maximum and minimum age bound in our dataset. For early postglacial and finiglacial eruptions the maximum age has to be tied to the inferred age of deglaciation of the area, based on available reconstructions of the ice sheet history (Geirsdóttir \BOthers., \APACyear2009; Patton \BOthers., \APACyear2017). The ages of subglacial eruptions are, in general, not as well constrained as those of the postglacial. Minimum age constraints for these eruptions are obtained from ice-sheet reconstructions and maximum ages are set to $30$ $\mathrm{ka}$. Helium-3 exposure ages and the geomorphological characteristics of the uppermost surface of tuyas can also be used to infer a chronology for a subset of subglacial eruptions, using the approach of Eason \BOthers. (\APACyear2015) as informed by the data of Licciardi \BOthers. (\APACyear2007).

A table of eruptions for which age, volume and chemical data is available is provided in the supplementary information. The information in this table is used to generate the plots provided for comparison with model results in this paper. The requirement of an unambiguous association between sample chemistry and eruption name, volume and age introduces some bias into our dataset: The lack of a clear link between the eruption name and chemistry means that our coverage of subglacial eruptions from the Reykjanes Peninsula is poor. Inevitably, erosion, superposition and lack of subsurface informaton introduce substantial uncertainties into any reconstruction of eruptive volumes.

We divide eruption units in WVZ further into WVZ-North (WVZN) and WVZ-South (WVZS) by latitude of $64^{\mathrm{o}}20^{\prime}0^{\prime\prime}$. Locations and types of eruption units of all the data we use here are plotted on the map in Figure 5.

The modelling-input distances of the four zones relative to the ice center shown in Table 3 are estimates with uncertainty of $\approx\pm 50$ $\mathrm{km}$ because the actual location of the ice center is unknown and also because of the uncertainty of the geometry of melt generation.

3.3 Eruption Types

In this section, we show how the melt ascent velocity can affect the relative volume proportion of different eruption types.

The modelling-input ice coverage radius as a function of time is known. We can therefore identify if an infinitesimal volume of melt that arrives at the surface at a particular location and time is erupted within the ice coverage radius or not. In other words, the model can divide eruptive volumes into subglacial group and subaerial group. The subglacial group corresponds approximately to the observational subglacial and finiglacial types combined. The subaerial group corresponds to the observational postglacial type. Our model does not divide the subglacial group further into subglacial and finiglacial types.

Figure 6 illustrates the model prediction that at faster melt transport (Figure 6b) there is a greater proportion of the subglacial volume (colored blue) compared to that at slower melt transport (Figure 6a). This is because faster melt transport will allow melts from depth to arrive at the surface sooner before the ice has gone. The sharp changes of subglacial to subaerial volume at $45$, $160$ and $180$ $\mathrm{km}$ are due to the three pauses of the glacial terminus at these three radial distances (Table 2).

The model also predicts that for the same time interval (such as $14.5$–[math] $\mathrm{kyrBP}$) the relative proportion of the subglacial volume to the total volume is dependent of the distance from the ice center. This is because while most of the melts in any location are produced over the same time period during deglaciation ($23.0$–$10.5$ $\mathrm{kyrBP}$), regions closer to the ice center remain covered by ice for a longer period of time. This allows a greater proportion of melts to arrive at the surface and erupt subglacially. The spatial dependence of the subglacial to subaerial volume ratio is also seen in observations. Figure 5 illustrates that in the regions closer to the center of Iceland there is a greater proportion of subglacial and finiglacial volumes (blue and green circles) than further out.

To compare our model results with the observations, we first note that the observational eruption volumes of units older than $14.5$ $\mathrm{kyrBP}$ are highly uncertain not only due to glacial erosion but also due to some older units are buried underneath younger eruptions. We therefore filter out eruptions older than $14.5$ $\mathrm{kyrBP}$ for both the model and the observational data. The model cumulative volumes at melt ascent velocity of $30$ and $1,000$ $\mathrm{m/yr}$ after the $14.5$ $\mathrm{kyrBP}$ filter are shown on Figure 6c and 6d. We use results from these two panels to calculate the model proportions of the subglacial volume ($=$ observational subglacial$+$finiglacial) and subaerial volume ($=$ observational postglacial) as shown on the two right bars of Figure 6e. For example, the bar plot of the model $30$ $\mathrm{m/yr}$ in NNVZ on Figure 6e has subglacial (blue) and subaerial (red) proportions equal to the subglacial (blue) and subaerial (red) plotting area proportions of Figure 6c in the x-axis range of $120$–$180$ $\mathrm{km}$.

We arrange the bar plot on Figure 6e from left to right by zone location from the closest to (WVZN) to the furthest from (REYK) the ice center. In each zone, the observational data has lower and upper estimates of subglacial volume due to age uncertainty of the subglacial units. The lower estimate (min. subglacial) shown on the left bar comes from the volume sum of the subglacial units with maximum age bound not exceeding $14.5$ $\mathrm{kyrBP}$. Whereas, the upper estimate (max. subglacial) comes from summing all the subglacial units with minimum age bound less than $14.5$ $\mathrm{kyrBP}$ (while the maximum age bound can exceed $14.5$ $\mathrm{kyrBP}$).

The model predictions for spatial dependence in the diachronous response agree well with the observational data. In REYK, the whole area is already ice-free by $14.5$ $\mathrm{kyrBP}$ and hence all the eruptions are subaerial. On the other hand, WVZN remains covered by ice over most of the time during the last deglaciation and so the majority of the eruption volumes are subglacial.

Results on Figure 6e also suggests that the melt ascent velocity is likely to be of the order of $100$ $\mathrm{m/yr}$. At below $30$ $\mathrm{m/yr}$, the subglacial volumes predicted by the model would be smaller than that of the observational lower bound estimates (min. subglacial). Nevertheless, we note that the model results depend on the distance along the ridge axis over which the melts are integrated (as estimated in Section 3.2). Similar to the model, the observational lava volumes in the four zones are also integrated over ridge lengths of $\sim 60$–$90$ $\mathrm{km}$.

3.4 Timing of the Peaks in Volcanic Productivity

Another way to estimate the melt ascent velocity is to use the timing of the peaks in volcanic productivity. On Figure 7, we plot the cumulative eruptive volume as a function of time for the model outputs and the observational data. This figure shows that the bursts in the cumulative lava volume predicted by the model at melt ascent velocities between $30$ and $1,000$ $\mathrm{m/yr}$ have timings approximately equal to that of the observations across all the volcanic zones to within the uncertainties of the lava ages and the modelling-input ice load history.

In the period during which the glacier terminus was sweeping through each zone (called transitional period), some areas in the zone are already ice-free while some areas are still covered by ice. This means that, in the transitional period, eruptions can be either subglacial or subaerial.

In the observational data sorted by age, the transitional period can be identified approximately by the period during which there are some alternations of the timeline orders between subglacial, finiglacial and postglacial types. The remaining two end periods are called subglacial and postglacial periods. The subglacial period consists of only subglacial type and the postglacial period consists of only postglacial type. In the model, the transitional period is identified by the period during which the ice radius is in the zone range (as listed in Table 3).

The timing of the periods in the model is controlled by the input ice history and the input zone range. Therefore, a well-matched timing of the periods between the model and the observational data helps us identify if the modelling-input of ice load and zone range reasonably reflect the actual values. The discrepancy between the observations and the model results at our preferred estimate of $100$ $\mathrm{m/yr}$ melt ascent velocity is likely to be because of the uncertainty of the observational eruption ages and the model axisymmetric ice assumption.

Nevertheless, the results in Figure 7 show that, at a melt ascent velocity of $30$ $\mathrm{m/yr}$ or below, the difference of the timing in the burst in volcanism between the model results and the observational data becomes significant. More quantitatively, with the average mantle melting depth of $\approx 50$ $\mathrm{km}$, reducing the melt ascent velocity from $100$ to $30$ $\mathrm{m/yr}$ will delay the burst timing by $\approx 1.2$ $\mathrm{kyr}$, which is significant compared to the uncertainty of the eruption ages. Our model implies a similar lower bound value to that of $50$ $\mathrm{m/yr}$ estimated in Maclennan \BOthers. (\APACyear2002) from the relative timing between the observed eruption ages and the deglaciation. An advantage of this work is that much broader geographical spread is accounted for. Our model examines eruptions in different volcanic zones; whereas, Maclennan \BOthers. (\APACyear2002) only examined eruptions in the NNVZ region.

3.5 Geochemical Response

3.5.1 Model Predictions

Figure 4b shows the volumetric supply rate of $\mathrm{La}$ to the crustal chambers normalized to the steady-state $\mathrm{La}$ concentration. Similarly to the whole melts in Figure 4a, the surge in the supply rate of $\mathrm{La}$ is delayed from the deglaciation period due to the finite speed of melt transport. However, while the volumetric melt supply rate curves (Figure 4a) are stretched in time, the $\mathrm{La}$ supply rate curves (Figure 4b) retain their shapes almost like the same time-series but time-shifted. This is because $\mathrm{La}$ is partitioned into melt at almost the same depth (near the solidus). Most of the $\mathrm{La}$ takes almost an equal time to arrive at the crust regardless of how fast the ascent rate is. Therefore, changing the ascent rate will not significantly spread the $\mathrm{La}$ flux out along the time axis. In contrast, melts are produced at different depths. They take different times to transport to the crust. The slower the ascent rate, the more the time delay between melts from different depths to arrive at the surface; hence, the more the spread of the melt supply rate curve along the time axis.

Figure 4c is the $\mathrm{La}$ concentration ($c_{\mathrm{La}}$) in the melt supply to the crustal chambers, which is equal to Figure 4b divided by Figure 4a. This would correspond to the concentration in the lava erupted at the surface if the magma mixing process in the crustal chamber were not present. The longer the magma is allowed to mix in the crustal chamber, the smaller the variation signal of the REE concentrations.

In the model, the effect of magma mixing on $c_{\mathrm{La}}$ in the lavas is equivalent mathematically to the time average concentration. The time period over which the average is performed is equal to the time duration that the magmas mix in the crustal chamber before they erupt. In a study of the chemical disequilibria between olivine, its melt-inclusions and the whole melts that surround the olivine in rock samples collected from Iceland, Maclennan (\APACyear2008) estimated that the magma residence time is of the order of a few hundreds to $\lesssim 1,000$ years. We therefore take the time average $c_{\mathrm{La}}$ over a period of $1,000$ years and the result is shown in Figure 4d. In other words, this panel is equal to the ratio between the $1000$-year standard moving average (SMA) of $\mathrm{La}$ volume supply to the crustal chamber ($1000$-year SMA of Figure 4b) and the $1000$-year SMA of the whole melt volume supply to the crustal chamber ($1000$-year SMA of Figure 4a).

Note that melt mixing also occurs "en-route" while melts are migrating from depths to the crustal chamber and meet together along the way. This geological process corresponds mathematically to the volume integration over the mantle melting domain $\mathcal{V}$ as shown in equation (9), taking into account the time delay due to the finite rate of melt transport $\Delta t$. In other words, the model results shown in Figures 4b, 4c and 4d as calculated by equations (9) and (10) have already taken into account the effect of en-route melt mixing.

Our results in Figure 4c and 4d show that the variation of $c_{\mathrm{La}}$ is strongly dependent on the melt ascent velocity. The lower the melt ascent velocity the higher the variation of $c_{\mathrm{La}}$. This effect can be explained as follows. During the deglaciation, the decompression rate in the mantle is maximum at the surface and decays exponentially with depth (as illustrated in Figure 3). This means that the extra melts generated during deglaciation are mostly produced at shallow depths in the mantle, which is depleted in $\mathrm{La}$. If the melt transport had been instantaneous, the extra melts produced at any depth at the same time would have travelled to the crust and mixed instantly and would have erupted at the surface with $\mathrm{La}$ depletion of up to $\approx 20\%$ as predicted by Jull \BBA McKenzie (\APACyear1996). In contrast, when the melt ascent velocity is finite, the extra melts produced at shallower depths during deglaciation will arrive at the surface before the extra melts produced at deeper depths. The slower the rate of melt transport the more likely the extra melts from shallow depths ($\mathrm{La}$ depleted) are to erupt before they mix with the extra melts from deep depths ($\mathrm{La}$ enriched). As a result, when the melt ascent velocity is sufficiently low, the first arrival of the extra melts produced during deglaciation will be much more depleted in $\mathrm{La}$ than that predicted by the instantaneous melt transport model of Jull \BBA McKenzie (\APACyear1996).

The eruptive $c_{\mathrm{La}}$ will recover back to near the steady-state concentration after the extra melts from the bottom of the melting region ($\mathrm{La}$ enriched) catch up and mix with the extra melts from shallow depths ($\mathrm{La}$ depleted) before they erupt. Moreover, the recovery of the eruptive $c_{\mathrm{La}}$ back to the steady-state will overshoot after the deglaciation ends. This phenomenon can be explained as follows. Once the deglaciation terminates, the glacially-induced decompression melting in the mantle will also terminate at all depths at the same time and the extra melt supply to the surface from shallow depths ($\mathrm{La}$ depleted) will run out before the extra melt supply from deep depths ($\mathrm{La}$ enriched). This is because the melts from greater depths take a longer time to arrive at the surface. Therefore, once the $\mathrm{La}$ depleted melt supply from the shallow depths runs out, the remaining majority of the erupted lavas will be the $\mathrm{La}$ enriched melts from deep depths and the eruption will become enriched in $\mathrm{La}$.

Figure 4d also shows that the timing of the periods during which the lavas are enriched or depleted in $\mathrm{La}$ is dependent on the melt ascent velocity. At slower melt ascent velocity, the peaks and the troughs of $c_{\mathrm{La}}$ are delayed further from the deglaciation periods. Hence, we can use this timing combined with the magnitude of the $c_{\mathrm{La}}$ variations to estimate the melt ascent velocity. We note that the $\mathrm{La}$ depleted lava volume dominates the $\mathrm{La}$ enriched lava volume. This can be seen in Figure 4. The troughs of $c_{\mathrm{La}}$ (Figure 4d) fall in the periods of the bursts in eruption rates (Figure 4a); whereas, the peaks of $c_{\mathrm{La}}$ (Figure 4d) fall outside those periods. Therefore, the $\mathrm{La}$-depletion signal is stronger than the $\mathrm{La}$-enrichment signal, which is also seen in observational data. The majority of the eruptions during the last deglaciation are depleted in $\mathrm{La}$.

3.5.2 Geological Observations

The eruptive $\mathrm{La}$ concentrations of observational data are from rock samples collected from Iceland by the previous studies (see Acknowledgements section and Supporting Information for details). These rock samples are of melts that have gone through fractionation/accumulation during cooling and crystallization processes in the crustal chambers. These processes modify the melt compositions from their original pre-crustal compositions.

We make fractionation/accumulation correction of $c_{\mathrm{La}}$ in each rock sample based on the $\mathrm{MgO}$ content of the sample. We assume that the pre-crustal melts have $14.0$ $\mathrm{wt\%}$ $\mathrm{MgO}$ as estimated in Maclennan \BOthers. (\APACyear2001). Rock samples that have $\mathrm{MgO}$ between $9.5$ and $14.0$ $\mathrm{wt\%}$ are assumed to have undergone crystallization of olivine-rich material with $40.0$ $\mathrm{wt\%}$ $\mathrm{MgO}$. If the melts underwent crystallization further below $9.5$ $\mathrm{wt\%}$ $\mathrm{MgO}$, they generate a gabbroic solid with $11.0$ $\mathrm{wt\%}$ $\mathrm{MgO}$. Some rock samples have higher $\mathrm{MgO}$ content than $14.0$ $\mathrm{wt\%}$ of the pre-crustal melts. We assume that these samples are from melts that have been influenced by the accumulation of olivine crystals with $40.0$ $\mathrm{wt\%}$ $\mathrm{MgO}$.

Due to the age uncertainty of eruption units, $c_{\mathrm{La}}$ cannot be plotted directly as that of the model in Figure 4. Most of the eruptions have their age estimated as a time band bounded by some geological events with identifiable age (e.g. tephra layers). In WVZ and NNVZ, most of the eruption units fall into one of the following age bands:

Glacial (pre $14.5$/$14.4$ $\mathrm{kyrBP}$) 2. 2.

Eruptive Pulse 1 ($14.5$/$14.4$ to $12.0$ $\mathrm{kyrBP}$) 3. 3.

Eruptive Pulse 2 ($12.0$ to $10.3$ $\mathrm{kyrBP}$) 4. 4.

Early Postglacial ($10.3$ to $8.9$/$8.0$ $\mathrm{kyrBP}$) 5. 5.

Steady-State Postglacial (post $8.9$/$8.0$ $\mathrm{kyrBP}$)

REYK zone is different in that it is the furthest from the ice center and became ice-free by $14.5$ $\mathrm{kyrBP}$, which is earlier than the other zones. Eruption units in REYK are divided into the following age bands:

Glacial (pre $14.5$ $\mathrm{kyrBP}$) 2. 2.

Early Postglacial 1 ($14.5$ to $13.0$ $\mathrm{kyrBP}$) 3. 3.

Early Postglacial 2 ($13.0$ to $10.2$ $\mathrm{kyrBP}$) 4. 4.

Steady-State Postglacial (post $10.2$ $\mathrm{kyrBP}$)

In each of these age bands, we take the volume-weighted average $\mathrm{La}$ concentration normalized to the Steady-State Postglacial $\mathrm{La}$ concentration and plot in Figure 8 for both the observational data and the model.

The right column of Figure 8 illustrates that different model melt ascent velocities result in different $\mathrm{La}$ concentration characteristics. In each age band, some values of melt ascent velocity may predict $\mathrm{La}$ depletion, whereas the others predict $\mathrm{La}$ enrichment. This is due to the effect of melt ascent velocity on the timing of the peaks and troughs of $c_{\mathrm{La}}$ (Figure 4d) that we discussed earlier.

The left column of Figure 8 shows that the model melt ascent velocity of $100$ $\mathrm{m/yr}$ yields similar $c_{\mathrm{La}}$ characteristics to that of the observations. At two extreme melt ascent velocities of $30$ and $1,000$ $\mathrm{m/yr}$, the $c_{\mathrm{La}}$ characteristics are significantly different from that of the observations both in the timing and the magnitude of the $c_{\mathrm{La}}$ variations.

In the Glacial age band, the observational $c_{\mathrm{La}}$ is elevated from the Steady-State Postglacial value in all the volcanic zones. This is likely to be due to the effect of glacial loading on depth-dependent melting suppression that occurred before the Last Glacial Maximum ($23$ $\mathrm{kyrBP}$). This feature is not included in our model here, which may explain why the model $c_{\mathrm{La}}$ in the Glacial age band is lower than that observed across all the volcanic zones. One of our future works will be to investigate this glacial loading effect.

3.6 Model Limitations

The accuracy of our results depends on several factors. The deviations of modelling input parameters from the actual geological values that are not well-constrained can be significant.

For example, the model $\mathrm{La}$ concentration is dependent on the time period over which the magma mixes in the crustal chamber. As mentioned in Section 3.5, the longer the magma residence time, the lower the variations of $\mathrm{La}$ concentrations. Also, the residence time may not be the same throughout Iceland as assumed in our model. A better constraint on the effective magma residence time in the chamber may therefore be required.

Our $~{}100$ $\mathrm{m/yr}$ estimate of the melt ascent velocity likely represents that of the melt produced during the GIA. At steady state, the decompression melting rate is significantly less. This leads to a significantly lower mass flux of melt and likely results in a slower rate of melt transport. This could be one reason why our ascent rate estimate is significantly higher than that in models of melt transport at Mid-Ocean Ridges (Burley \BBA Katz, \APACyear2015; Crowley \BOthers., \APACyear2015). A melt ascent velocity of $100$ $\mathrm{m/yr}$ is $\sim 2$ orders of magnitude faster than would be predicted from simple models of diffuse porous flow. As has been noted in several previous studies, such rapid melt ascent velocities suggest that some focusing or channelization of melt must occur during transport (e.g. Kelemen \BOthers. (\APACyear1997)).

In more elaborate fluid dynamic models (e.g. McKenzie (\APACyear1984)), the melt velocity varies with depth. Melt flow starts slow at the base of the melting region and ascends at a faster rate as it migrates to a shallower depth where the porosity is higher. Therefore, in any region below the GIA average melting depth, the melt ascent velocity is likely to be below our estimate. This depth-dependent melt ascent velocity would cause a more time delay of the burst of the $\mathrm{La}$ supply rate to the crustal chamber than that predicted by our model in Figure 4b. A larger time delay between the melt supply and the $\mathrm{La}$ supply would increase the time intervals during which the $\mathrm{La}$ concentration (Figures 4c and 4d) is depleted or enriched.

2-D fluid dynamic models of melt and trace element transport (e.g. Spiegelman (\APACyear1996)) also predict an across-axis variation in the erupted melt composition. In our model, we assume complete melt mixing and extraction on the ridge axis. This produces only a single average concentration of $\mathrm{La}$ at each snapshot in time. Spiegelman (\APACyear1996) also showed that the convergence of melt to the ridge axis in passive ridge flow leads to an enrichment of incompatible elements in the erupted melt by almost a factor of $2$ (for $D^{i}\leq 0.01$) from that in the 1-D column model. If the full solution of melt transport had been incorporated into our model, the $\mathrm{La}$ concentration at steady state would have also been increased by a factor of $\sim 2$. If the same enrichment factor ($\sim 2$) also uniformly applies to that during the GIA, our model results of the $\mathrm{La}$ concentration (normalized to the steady state value) would remain unchanged. However, the flow fields of the solid mantle and of the melts during the GIA are certainly different from those at steady state. Therefore, the enrichment factor during the GIA is unlikely to be uniformly the same as that at the steady state. How much the enrichment factor varies still remains to be explored. In our future work, we would like to incorporate full melt transport solutions into the model to understand how good the constant ascent rate approximation is.

The real ice sheet shape may be significantly deviated from the axisymmetric shape that we use. While our axisymmetric assumption helps simplify the computations, a modelling-input ice sheet with more detailed 3D shape may have an important role in controlling the accuracy of the model results.

Last but not least, the time evolution of the ice sheet shape we input into our model may be significantly different from the actual ice sheet. For example, the glacier may extensively re-advance in some periods during the last deglaciation. Glacial advance will increase the load on the surface, which will lead to pressure increase in the mantle. This will suppress mantle melting and can also affect the REE concentrations as discussed in Section 3.5. The modelling of the effects of glacial advance on mantle melting beneath Iceland may therefore be important.

4 Conclusions

The consequences of a finite melt ascent velocity on lavas erupted during the last deglaciation are:

Volume proportions of different eruption types: Faster melt transport will allow more melts to arrive at the surface and erupt sooner when the ice is still present. This means that there will be a greater proportion of subglacial and finiglacial volumes relative to postglacial volume. 2. 2.

Relative timing between the bursts in the eruption rates and the deglaciation: Higher melt ascent velocity will transport the extra melts produced during deglaciation to the surface faster. This will result in a smaller time-lag between the bursts in the eruption rates and the deglaciation. 3. 3.

Variations of REE concentrations: Slower melt ascent velocity will result in a greater time-lag between melts from shallow depth (REE depleted) and melts from deep depth (REE enriched) arriving at the surface. This will cause higher variations of REE concentrations in the lavas.

Our numerical model estimates that the Icelandic melt transport from the upper mantle melting region to the surface during the last-deglaciation has an average melt ascent velocity of the order of $\sim 100$ $\mathrm{m/yr}$.

Appendix A Mantle Flow

In steady state, the velocity components of the mantle flow follow the corner flow solutions. In Cartesian coordinates, they can be written as (Batchelor (\APACyear2000); Spiegelman \BBA McKenzie (\APACyear1987))

[TABLE]

where

[TABLE]

$v_{x}$ is the horizontal velocity component perpendicular to the ridge, $v_{y}$ is the horizontal velocity component parallel to the ridge, $v_{z}$ is the velocity component in the vertical direction, $U_{0}$ is the half-spreading velocity of the ridge and $\alpha$ is the ridge angle from the horizontal.

Ice load change causes glacially-induced isostatic adjustment (GIA). In cylindrical coordinates, the semi-analytical solutions to the axi-symmetric GIA response in the viscous half-space mantle are

[TABLE]

where

[TABLE]

$v_{r}$ is the radial component of the velocity, $v_{\theta}$ is the azimuthal component of the velocity, $v_{z}$ is the vertical component of the velocity, $P$ is the pressure in the mantle, $\rho_{i}$ is the density of ice, $\rho_{s}$ is the density of the mantle, $\dot{h}$ is the time-derivative of the thickness of ice sheet, $\mathcal{H}_{n}[f]$ is the $n^{\text{th}}$-order Hankel transform of function $f$, $\mathcal{H}^{-1}_{n}[F]$ is the $n^{\text{th}}$-order inverse Hankel transform of function $F$, $\mathrm{J}_{n}$ is the $n^{\text{th}}$-order Bessel function of the first kind and $k$ is the wavenumber. $t_{\mathrm{iso}}$ in the $\widetilde{w}\,(k,t)$ expression is the time at which the mantle is in isostatic equilibrium.

Equation (2) shows that the glacially induced decompression rate depends on the history of the deglaciation rate $\dot{h}$ and is attenuated exponentially with depth by the $e^{kz}$ factor ($z<0$ in the mantle). For the ice sheet shape in equation (1), the Hankel transform of the rate of change of ice load is analytical:

[TABLE]

Appendix B Numerical Methods

Calculations of the melting rates, the eruption rates and the REE concentrations (equations (5), (6), (8) and (9)) require temporal and spatial integrations over a finite domain. We perform these numerical integrations using the trapezoidal rule. We discretize the spatial domain using a 3-D rectangular grid with uniform horizontal and vertical resolutions of $5$-by-$5$ $\mathrm{km^{2}}$ and $0.5$ $\mathrm{km}$ respectively.

For an REE with $D^{i}\ll 1$, the concentration $c^{i}_{l}$ changes rapidly with depth near the solidus. As can be seen in equation (7), $c^{i}_{l}$ drops sharply with $F$ near the solidus $F=0$. By adopting the trapezoidal rule to integrate equation (8) along the depth using values of $c^{i}_{l}$ at the grid vertices alone, the trapezoidal error can be very significant. In order to resolve this rapid change within each cell of the grid, we calculate $c^{i}_{l}$ inside the cell using equation (7) with $F=F(z)$ that is obtained from linear interpolation of the face-averaged $F$ between the lower face and upper face of the cell with depth $z$. The face-averaged value of a face is simply the mean value of the four vertices at the corners of the face. This technique helps improve the model $c^{i}_{l}$ accuracy significantly.

Calculations of mantle flow and decompression rates such as equation (2) involve the inverse Hankel transform, which requires numerical integration of the wavenumber $k$ from [math] to $\infty$. All the mathematical expressions in the model that require inverse Hankel transform contain an attenuation factor $e^{kz}$, which decays exponentially with the wavenumber $k$ in the mantle ($z<0$). Therefore, the numerical integration of the inverse Hankel transform from $k=0$ to $\infty$ can be truncated when the attenuation factor $e^{kz}$ is negligibly small. In our model, the melting region is at $z=-20$ $\mathrm{km}$ or below. We truncate the integration at $k=2/3$ $\mathrm{km}^{-1}$, which corresponds to $e^{kz}\sim 2\times 10^{-6}$ or below.

The variations of integrands of all the inverse Hankel transform involved in the model are dominated by the ice load function in the $k$-domain (equation (3)). This function consists of the $n^{\text{th}}$-order Bessel functions of the first kind $\mathrm{J}_{n}$ with $n=\pm 1/4$ and $\pm 3/4$, all of which have the same argument $=kr_{m}/2$ where $r_{m}$ is the ice radius. This means that the ice load function in the $k$-domain varies with $k$ at a frequency of $\sim r_{m}/2\sim 100$ $\mathrm{km}$. We therefore use the trapezoidal strip size $\mathrm{d}k=1/1440$ $\mathrm{km}^{-1}$, which gives $\mathrm{d}kr_{m}/2\leq 10^{-1}$. This corresponds to having at least $N_{k}=10$ trapezoidal strips per unit length in the non-dimensional $k$-domain since $N_{k}=2/(\mathrm{d}kr_{m})\geq 10$.

The decompression rate $=-DP/Dt$ at time $t$ can be calculated directly from equation (2) independently from any information in the previous time steps. This means that the time-step size does not affect the accuracy of the model. We use a uniform time-step size of $50$ $\mathrm{yr}$.

Acknowledgements

We would like to thank Marc W. Spiegelman, Deborah Eason and an anonymous reviewer for their constructive and thorough reviews, which greatly helped improve this manuscript. This research is funded by the Cambridge Trust and the Leverhulme Trust. All the data of rock samples we use (as provided in the Supporting Information) come from the following published sources: Brandon \BOthers. (\APACyear2000); Brandon \BOthers. (\APACyear2007); Breddam \BOthers. (\APACyear2000); Burnard \BBA Harrison (\APACyear2005); Chauvel \BBA Hémond (\APACyear2000); Condomines \BOthers. (\APACyear1983); Debaille \BOthers. (\APACyear2009); Dixon \BOthers. (\APACyear2000); Dixon (\APACyear2003); Eason \BBA Sinton (\APACyear2009); Eason \BOthers. (\APACyear2015); Elliott \BOthers. (\APACyear1991); Fitton \BOthers. (\APACyear2003); Füri \BOthers. (\APACyear2010); Gee, Taylor\BCBL \BOthers. (\APACyear1998); Gee, Thirlwall\BCBL \BOthers. (\APACyear1998); Geirsdóttir \BOthers. (\APACyear2009); Hardarson \BOthers. (\APACyear1997); Hemond \BOthers. (\APACyear1993); Jakobsson \BOthers. (\APACyear1978); Jónasson (\APACyear2005); Kempton \BOthers. (\APACyear2000); Kokfelt \BOthers. (\APACyear2006); Koornneef \BOthers. (\APACyear2012); Kurz \BOthers. (\APACyear1985); Maclennan \BOthers. (\APACyear2001); Maclennan \BOthers. (\APACyear2003); Maclennan \BOthers. (\APACyear2004); Macpherson \BOthers. (\APACyear2005); Nicholson \BOthers. (\APACyear1991); Nielsen \BOthers. (\APACyear2007); Peate \BOthers. (\APACyear2009); Peate \BOthers. (\APACyear2010); Poreda \BOthers. (\APACyear1986); Sæmundsson (\APACyear1991); Sæmundsson \BOthers. (\APACyear2012); Sæmundsson \BOthers. (\APACyear2016); Sigurdsson \BOthers. (\APACyear1978); Sinton \BOthers. (\APACyear2005); Skovgaard \BOthers. (\APACyear2001); Slater \BOthers. (\APACyear2001); Sobolev \BOthers. (\APACyear2008); Sólnes \BOthers. (\APACyear2013); Stracke \BOthers. (\APACyear2003); Thirlwall \BOthers. (\APACyear2004); Thirlwall \BOthers. (\APACyear2006)