Thermal impact of magmatism in subduction zones

TL;DR

This paper demonstrates that magma transport significantly influences the thermal structure of subduction zones, explaining observed heat flow and temperature discrepancies with traditional models, and affecting volcanic activity.

Contribution

It introduces numerical models showing magma transport's role in altering subduction zone thermal structures, resolving previous temperature prediction discrepancies.

Findings

Magma transport raises local temperatures by ~300 K.

Heat transported by magma explains higher observed heat flow.

Magma advection dominates latent heat deposition.

Abstract

Magmatism in subduction zones builds continental crust and causes most of Earth's subaerial volcanism. The production rate and composition of magmas are controlled by the thermal structure of subduction zones. A range of geochemical and heat flow evidence has recently converged to indicate that subduction zones are hotter at lithospheric depths beneath the arc than predicted by canonical thermomechanical models, which neglect magmatism. We show that this discrepancy can be resolved by consideration of the heat transported by magma. In our one- and two-dimensional numerical models and scaling analysis, magmatic transport of sensible and latent heat locally alters the thermal structure of canonical models by $\sim$300 K, increasing predicted surface heat flow and mid-lithospheric temperatures to observed values. We find the advection of sensible heat to be larger than the deposition of latent heat. Based on these results we conclude that thermal transport by magma migration affects the chemistry and the location of arc volcanoes.

Equations

$$\frac{\partial T}{\partial t}+{\mathbf{v}}_s\cdot \nabla T+{\mathbf{v}}_D\cdot \nabla T=\kappa {\nabla }^{2}T-\frac{L}{\rho c_p}\Gamma ,$$

$$\rho c_p\overline{W_0}\frac{dT}{dz}-\rho c_p{\Psi }^{\ast }=\frac{d}{dz}\left(\rho c_p\kappa \frac{dT}{dz}\right)-L\Gamma ,$$

$$\text{Pe}{T}^{\prime }-\Psi =T^{\prime \prime }-\text{Pe}\text{St}(T^{\prime }+\Delta T_H),$$

$$Q\approx \frac{F_V\rho c_p\Delta T}{A}\approx 80{mW/m}^2,$$

$$R\approx \frac{\rho c_p\left|{\mathbf{v}}_D\right|(\Delta T/H)}{L\Gamma }\approx \frac{c_p\Delta T}{L}\frac{\rho \left|{\mathbf{v}}_D\right|}{\Gamma H}\approx \frac{c_p\Delta T}{L}\approx 3.2,$$

$$T_s=T_{s0}-\rho gz/\gamma -M_2{c}_2^s-M_3{c}_3^s.$$

$$c_2^l$$

$$c_3^l$$

$$F=\frac{\widetilde{c_{0j}}-c_j^s}{c_j^l-c_j^s}.$$

$$F=\frac{\widetilde{c_{03}}-c_3^s}{\Delta c_3}=\frac{T-T_{s0}+\rho gz/\gamma +M_2\widetilde{c_{02}}+M_3\widetilde{c_{03}}}{M_2\Delta c_2+M_3\Delta c_3}.$$

$$\frac{\partial F}{\partial T}=\frac{1}{M_2\Delta c_2+M_3\Delta c_3}.$$

$${\frac{\partial F}{\partial T}|}_{F=0}<{\frac{\partial F}{\partial T}|}_{F=1}.$$

$$c_3^s=D{c}_3^l\Rightarrow \Delta c_3=c_3^s(1/D-1),$$

$$T-$$

$$=(M_2\Delta c_2/\Delta c_3+M_3)(\widetilde{c_{03}}-c_3^s),$$

$${\frac{\partial F}{\partial T}|}_{F=0}$$

$${\frac{\partial F}{\partial T}|}_{F=1}$$

$$M_3\widetilde{c_{03}}D(1-D)\ll M_2\Delta c_2\ll M_3\widetilde{c_{03}}(1/D-1).$$

$$T(F=0)$$

$$T(F=1)$$

$$X_{{\mathrm{H}}_{2}0}^{\mathrm{sat}}=12.00P^{0.6}+1.00P,$$

$$T_s=T_{s0}-\rho gz/\gamma -M_2{c}_2^s-M_4(c_3^s-c_{\mathrm{sat}}^s)-M_3{c}_{\mathrm{sat}}^s,$$

$$\frac{d}{dz}(\varphi w^l)$$

$$\frac{d}{dz}\overline{w}$$

$$\frac{\varphi w^l}{\overline{W_0}}+\frac{(1-\varphi )w^s}{\overline{W_0}}=1,$$

$$\frac{d}{dz}(\varphi w^l{c}_j^l)$$

$$\frac{d}{dz}\overline{w{c}_j}$$

$$F{c}_j^l+(1-F)c_j^s=\frac{\overline{W_0{c}_{0j}}}{\overline{W_0}}\equiv \widetilde{c_{0j}}.$$

$$\rho c_p\overline{W_0}\frac{dT}{dz}=-L\Gamma +\frac{d}{dz}\left(\rho c_p\kappa \frac{dT}{dz}\right)+\rho c_p\Psi ,$$

$$T^{\prime \prime }=-\Psi +\text{Pe}\left[T^{\prime }(1+\text{St})+\text{St}\Delta T_H\right],(0\le z\le 1).$$

$$\tilde{\Gamma }\equiv \frac{H^2}{\kappa }\frac{\Gamma }{\rho }=\text{Pe}\frac{\partial F}{\partial T}(T^{\prime }+\Delta T_H),$$

Full text

I Introduction

Petrological estimates of sub-arc temperature conditions in both continental and oceanic subduction zones are systematically higher than predicted by thermal models, typically by 200–300 K, at depths less than $\sim$70 km Kelemen2003 ; Perrin2016 . Similarly, measurements of geothermal heat flow in SW Oregon and NE Japan are higher than predicted by approximately 50–100 mW/m2 near the volcanic arc Kelemen2003 ; Furukawa1993 . Geophysical evidence from seismic and magnetotelluric imaging of high temperatures and/or magma at depth under volcanic arcs Zhao2007 ; Syracuse2008 ; Rychert2008 ; McGary2014 is consistent with the emerging consensus that the shallow arc temperatures in subduction zones are hotter than canonical models predict.

In canonical models, the thermal structure of subduction zones is calculated as a balance between thermal diffusion and advection. Heat is advected by the creeping solid mantle flow within the wedge-shaped region between the subducting slab and overriding lithosphere McKenzie1969 . Previous modelling efforts to resolve the discrepancy with observations have involved varying the prescribed geometry of subduction, the coupling between mantle and slab, and the rheological model of the mantle Kelemen2003 ; Furukawa1993 . Inclusion of frictional heating along the slab top in the seismogenic zone increases heat flow in the fore-arc Gao2014 . None of these efforts have been successful in explaining both the amplitude of the thermal observations and their position relative to the volcanic arc.

It is known that hydrous fluids are released from the subducting slab by de-volatilization reactions Schmidt2014 and percolate upward into the mantle wedge. There they reduce the solidus temperature, promote melting, and hence become silicic as they ascend. During their ascent, the magmas traverse from cooler mantle adjacent to the slab, to hotter mantle at the core of the wedge, to cooler mantle at the base of the lithosphere. They advect heat between these regions and consume or supply latent heat with melting and freezing. Despite the copious production of magma in subduction zones, these thermal processes have been neglected from almost all previous models. One exception, a scaling argument comparing advective heat transport by magma flow to thermal diffusion, suggests that magma flow may be significant Peacock1990 . Similarly, hydrothermal circulation in the crust may play a role in cooling the slab in the fore-arc region Spinelli2016 . In this paper we assess the role of magmatic processes in altering the thermal structure of the wedge and lithosphere. Our approach is based on theory for two-phase dynamics of the magma–mantle system McKenzie1984 . We quantify the magmatic transport of sensible and latent heat, focusing on the physical mechanisms and their controls, rather than on any particular subduction zone.

II Methodology

Magma migration in the mantle is a two-phase flow, governed by continuum equations of mass and momentum conservation for the solid (mantle) and melt (magma) McKenzie1984 ; Rudge2011 . The thermal and compositional structure is governed by equations of conservation of energy and chemical species. Our approach is to prescribe the magmatic flux and investigate how the thermal structure responds. This response is determined from energy conservation in the form of a heat equation:

[TABLE]

$T$ denotes temperature, $t$ time, $\kappa$ thermal diffusivity, $\rho$ density, $c_{p}$ specific heat capacity, $L$ latent heat, and $\Gamma$ melting rate. We neglect differences between the thermal properties of the phases because these do not affect the solution at leading order. The velocity variables involved are: solid mantle velocity $\boldsymbol{v}_{s}$, liquid magma velocity $\boldsymbol{v}_{l}$, the Darcy (or segregation) flux $\boldsymbol{v}_{D}\equiv\phi(\boldsymbol{v}_{l}-\boldsymbol{v}_{s})$, where $\phi$ is the porosity.

In the absence of magma, $\boldsymbol{v}_{D}=0$ and $\Gamma=0$ and eqn. (1) reduces to the heat equation used in canonical mantle convection calculations. In the presence of magma, two relevant terms are non-zero: first, an advective term associated with the segregation flux of magma $\boldsymbol{v}_{D}$; second, a latent heat sink associated with melting ($\Gamma>0$), which becomes a source in the case of freezing ($\Gamma<0$). The petrological model for $\Gamma$ is described in Sec. S1, Supplementary Material, and was inspired by previous studies of mantle melting in the presence of water Hirschmann1999III ; Katz2003 ; Keller2016 .

By the considerations above and the results below, we emphasize that the latent heat of phase change is not the only thermal contribution from magmatism; there is also advective transport by the magma. In what follows, we consider the relative importance of these mechanisms.

III Results

III.1 One-dimensional model

So-called ‘melting-column models’ have been used to understand mid-ocean ridge magmatism, where the main cause of melting is decompression of the upwelling mantle Ribe1985 ; Asimow1999 ; Hewitt2010 . Subduction zones are a considerably more complex environment, but we adapt ideas from melting-column models to investigate how magmatism modifies their thermal structure. The column model is fully derived and described in more detail in Sec. S2, Supplementary Material. A one-dimensional, steady-state heat equation can be written

[TABLE]

where $\Psi^{*}$ is the dimensional version of the source term, discussed below. We rescale lengths by the height of the column $H$, velocities by the diffusive scale $\kappa/H$, and $\Psi^{*}$ by $\kappa/H^{2}$. Then eqn. (2) becomes

[TABLE]

where $\Psi$ is the rescaled version of the source term, discussed below. $\Delta T_{H}$ is the adiabatic temperature drop between slab and surface; primes denote a derivative with respect to position (e.g., $T^{\prime}$ is a rescaled vertical temperature gradient). Two dimensionless numbers control the behaviour of the system: a Péclet number $\textrm{Pe}\,=H\overline{W_{0}}/\kappa$ is the scaled volume flux at the base of the column; a Stefan number $\textrm{St}\,=(L/c_{p})\partial F/\partial T$ is the scaled isobaric productivity that quantifies the ratio of latent to sensible heat ($F$ is the degree of melting). Hydrous flux melting has low isobaric productivity Hirschmann1999III so the Stefan number is small.

The mantle flow in subduction zones is far from one-dimensional; a corner flow is driven by the motion of the subducting slab McKenzie1969 . A key step in representing corner flow in a column model is to introduce a spatially variable, volumetric heating term $\Psi$ that mimics the effects of large-scale mantle flow, which tends to supply heat into the column. We infer $\Psi$ from a single-phase, two-dimensional thermomechanical reference model that is shown in Figure 1; the domain geometry and temperature-dependence of viscosity are as given in a study that outlined broadly representative models of subduction vanKeken2008 . From the reference model, we extract a vertical temperature profile at some position of interest $T_{\mathrm{ref.}}(z)$ and use it to calculate the source term $\Psi=-T^{\prime\prime}_{\mathrm{ref.}}$. The source term is constructed such that the solution of equation (3) in the absence of magma flow ($\textrm{Pe}\,=0$) is $T=T_{\textrm{ref.}}$, i.e., the single-phase result. For $\textrm{Pe}\,>0$, this approach is reasonable provided melt does not drastically change the large-scale mantle dynamics, a prospect we consider later.

Figure 2 shows results of the 1D column calculations. These are obtained for the column rising from slab where it is 100 km deep. This choice is roughly consistent with the observed mean depth beneath arc volcanoes EnglandEngdahl2004 ; Syracuse2006 . The flux at the base of the column is varied within the range suggested by a previous study Wilson2014 . Dimensionally, this range corresponds to fluxes between 0.2–2 m/kyr. Panel (a) shows profiles of the absolute temperature; panel (b) shows the temperature difference compared to the single-phase (magma-free) reference case. The change in temperature from the reference state increases with the imposed flux and is significant even at the lower end of the plausible range Wilson2014 . Immediately above the slab, upward flow reduces the mantle temperature as material is transported from the relatively cold slab. Nearer the surface, the effect is reversed as upward flow brings warm material from the mantle into the lithosphere. This effect is supplemented by latent heat associated with melting and solidification, shown in panel (c). Above the slab, melting of the mantle wedge facilitated by the presence of water consumes latent heat. Nearer the surface, solidification of the melt deposits latent heat. The maximum degree of melting (d) is increased because of the elevated temperatures, which will have a significant geochemical signature Turner2016 . It is interesting to note that the maximum degree of melting does not vary monotonically, but peaks at an intermediate Péclet number between 2 and 5.

The main physical mechanism giving rise to this thermal response is advection by the magma; latent heat release reinforces the advective heat flux. Additional calculations, shown in Figure 3, demonstrate that latent heat has a subordinate effect on the temperature profiles. Other calculations shown in Fig. 3 indicate that these results are robust to changes in the parameterization of hydrous flux melting (either to mimic more closely a more detailed parameterization Katz2003 , or by arbitrarily doubling the Stefan number). The relative importance of latent to specific heat is controlled by the Stefan number St . This is typically relatively small; $\textrm{St}\,<0.1$ throughout the temperature range encountered ($T_{\mathrm{ref.}}\leq 1250^{\circ}$C, above a slab 100 km deep), as shown in Supplementary Material, Fig. S3. If the Stefan number were much larger, latent heat release would be comparable to thermal advection by magma (Fig. 3). A larger Stefan number may be relevant for magmatic environments dominated by melting at high isentropic productivity, above the anhydrous solidus (i.e., plumes and mid-ocean ridges). But subduction zones are characterized by low-productivity hydrous-flux melting Hirschmann1999III , associated with a small Stefan number, and hence the role of latent heat is relatively minor.

III.2 Two-dimensional thermal model with magma migration

Two-dimensional effects that are neglected in column models, such as lateral diffusion and changes to viscosity structure and mantle flow, require a more careful treatment. We next consider the thermal consequences of magmatic advection by modifying a canonical, two-dimensional reference simulation of a subduction zone vanKeken2008 to include a prescribed segregation flux $\boldsymbol{v}_{D}$ in the heat eqn. (1). We assume that magma segregates purely vertically, driven by the density difference between solid and liquid phases. We prescribe this flow in terms of Gaussian profiles centred at the typical position of the arc volcano EnglandEngdahl2004 ; Syracuse2006 . Our numerical scheme solves iteratively for thermal structure and solid flow, which are fully coupled through advection and the temperature dependence of mantle viscosity, until a steady state is achieved. The thermal impact of magmatism is then defined as the difference between the calculated and reference temperature fields.

The two-dimensional calculations, shown in Figure 4, predict that magmatic transport substantially alters the thermal structure in subduction zones. The main effect is to raise temperatures near the base of the lithosphere, where warm material is transported from the mantle upward. These 2D results are qualitatively consistent with the 1D column models (cooling above the slab, warming near the surface), indicating that the physical mechanisms discussed above remain pertinent. However, some features only occur in two dimensions. For example, cooling is observed immediately above the slab-top deeper than 100 km; this is caused by advection with the mantle flow. Thus the thermal impact of magmatism is distributed beyond the imposed region where the magma flows.

Our standard estimate of the magmatic flux uses a Gaussian velocity profile (Fig. 4d) with a peak velocity of 2 m/kyr and a width of 10 km, giving a total flux comparable to global estimates Reymer1984 ; Crisp1984 ; England2010 . In this case, magmatism raises temperatures by up to 270 K (Fig. 4b). We also consider a magmatic flux 50% smaller or larger than this standard case. Temperatures are raised by $\sim$150 K (Fig. 4a) with the lower estimate. The higher estimate raises temperatures by up to 380 K (Fig. 4c). In three dimensions, the thermal effect local to arc volcanoes would likely be even greater due to along-strike flow focussing.

Figure 5 shows the results of additional calculations that explore the sensitivity to different parameter values and modelling choices that are consistent with observational constraints. For all these calculations, we compare against the standard magma flux case (Fig. 4b). For the impatient: these sensitivity experiments show that our key conclusion — that magmatism has a significant thermal effect — is robust.

First, we find that the total magma flux is more significant that the width of the flow. In Model Experiment 1, we show that similar temperatures are obtained by doubling of the width of the magma flow while halving of its magnitude to hold the total flux constant. The wider flow has a slightly lower peak (by 40 K) and is slightly more diffuse. However, these differences are minor compared to those associated with varying the total magma flux (Fig. 4a,c). The width of the thermal response is controlled primarily by the balance between advective heat transport by the magma and thermal diffusion.

Second, we consider the effect of the viscous coupling between the solid velocity and the temperature field. We partially decouple the model by holding the solid velocity field fixed at the reference conditions associated with the reference temperature field (i.e., that without magmatism). In Experiment 2, we show that the semi-decoupled calculations have a significantly smaller thermal response. The mechanism is as follows: in the fully coupled calculations, the elevated temperatures caused by magmatism lower the mantle viscosity, increasing the mantle wedge circulation, which is shown in Figure 6. This leads to increased heat transport toward the arc (a positive feedback). The effect of coupling is more pronounced with smaller plate thickness because there is a larger region of mantle flow where the viscosity is reduced, leading to faster circulation (cf. Exps. 3a and 3b in Fig. 5).

Third, we consider the effect of the imposed thickness of the overiding plate (Exp. 3a of Fig. 5). The thermal effect of magmatism decreases slightly with increasing plate thickness. This is associated with cooler temperatures in the reference state, reducing the advection of heat by the magma. The decrease is also aided by the fact that the coupling to the solid velocity becomes a less significant positive feedback as plate thickness increases (Exp. 3b of Fig. 5, which is relatively similar to Exp. 3a).

Fourth, we consider the effect of slab–wedge coupling (Exp. 4 of Fig. 5). We increase the slab–wedge coupling depth from 50 km to 80 km, a value suggested by Ref. Wada2009 on the basis of fore-arc heat flow measurements. This has a significant effect on the reference state without magmatism. However, it has only a small effect on the thermal effect of magmatism itself.

Fifth, we consider the effect of slab dip (Exp. 5 of Fig. 5). We double the slab slope from 1:1 to 2:1. Again, we find that the thermal effect of magmatism is qualitatively very similar to the standard case in Fig. 4b.

Finally, in Figure 7, we consider the transient evolution towards steady state. We use an initial condition corresponding to old oceanic lithosphere and impose the same fixed magma flux. The thermal effect of magmatism evolves to a steady state over a period of about 50 Myr, controlled by thermal diffusion, although the thermal structure much further away from the arc evolves on a longer timescale Hall2012 . The transient spatial pattern of elevated sub-arc temperatures is consistent with the steady-state pattern. However, the magnitude of the thermal effect depends on the age of the subduction zone.

In the Supplementary Material, Sec. S3, we consider separately the magmatism associated with each of the major slab dehydration reactions that occur at various depths.

In summary, in each sensitivity test, we find that although small quantitative differences in the results are produced, the overall behaviour and the basic conclusion is similar. Thus the thermal effect of magmatism we show in Figure 4 is robust; the details will vary between subduction zones, but the physical effect is to significantly modify the thermal structure from that predicted by canonical models.

IV Discussion and Conclusions

Our results are consistent with heat flow and petrological observations. The elevated heat flow measured in subduction zones, shown in Figure 8, can be associated with elevated near-surface temperatures. This elevated heat flow is strongest at the position of the arc, over a width of around 50 km. The width is determined by thermal diffusion rather than the imposed width of magma flow. Our models that use a magma flux between the standard and high values are consistent with heat-flow observations near the volcanic arc. Note that the low fore-arc heat flow in our models is an artefact of the simplified geometry, particularly the constant slab dip. Furthermore, hydrothermal circulation in the subducting crust has a significant thermal effect in the fore-arc region, consistent with heat flow observations along the Chilean subduction zone Spinelli2016 . Similarly, we find that magmatic flow has a significant thermal effect in the sub-arc region, consistent with heat flow observations there.

Evidence from petrological observations in Figure 9 suggests that temperatures in subduction zones are some 200–300 K hotter than would be expected on the basis of canonical models of mantle flow alone Kelemen2003 ; Perrin2016 . This discrepancy peaks at around 60 km depth, comparable to the depth where we find magmatism has the greatest thermal impact. Inclusion of melt migration in thermal models can reconcile much of this discrepancy. This consistency between observation and thermal modelling supports the hypothesis that magmatism significantly alters the thermal structure of subduction zones.

Scaling arguments also support our hypothesis. Indeed, it is possible to approximate the effect on heat flow due to magmatic advection as follows. The elevated heat flow is

[TABLE]

based on a global magma flux $F_{V}=1$ km3 yr*-1* Reymer1984 , density $\rho=3\times 10^{3}$ kg m*-3*, heat capacity $c_{p}=1.2\times 10^{3}$ J kg*-1* K*-1*, $\Delta T\approx 1350$ K, and an area of elevated heat flow $A\sim 2\times 10^{12}$ m2 (the total length of $50\times 10^{3}$ km and an assumed width of 40 km). This is consistent with Fig. 8.

We can also estimate the ratio $R$ of advective heat transport by magma to the latent heat release (the two mechanisms by which magmatism changes the thermal structure):

[TABLE]

where $L=5\times 10^{5}$ J kg*-1*. We used the fact that ${\rho\left|\boldsymbol{v}_{D}\right|}/{\Gamma H}\approx 1$ on average at steady state, since there is a balance between melt production, melt extraction, and melt solidification. Therefore, magmatism has a significant thermal effect and this effect is mainly due to advection by the magma. This latter finding is in contrast to a previous, simpler, one-dimensional model England2010 ; Perrin2016 .

The thermal signature of melt migration should be considered when interpreting heat flow, petrologic, gravity, and seismic data. Seismic velocities and attenuation depend strongly on temperature Takei2017 . Thus our results suggest that a part of the measured low seismic velocities and high attenuation beneath the arc is likely associated with high temperatures. However, the relatively small spatial extent of the thermal anomalies we predict ($\sim$50–100 km) will make them difficult to observe seismically. A perturbation as large as 300 K also increases the maximum degree of melting, which in turn affects the chemistry of arc volcanoes (or our inferences about the mantle made on the basis of geochemical measurements) Turner2016 . It also significantly affects the solid mantle flow through reduction of mantle viscosity, leading to increased circulation in the mantle wedge vanKeken2002 . Furthermore, thermal structure affects magma pathways in subduction zones, focussing magmas along the thermal lithosphere from a broader area to beneath the arc volcanoes Sparks1991 ; Wilson2014 . Thus, coupled mantle–magma flow may well affect the location of arc volcanoes themselves, consistent with evidence from global systematics England2010 .

Author contributions

R.F.K. conceived the study. D.R.J. and M.T. developed the one-dimensional melting column model and petrological model of melting. D.R.J., R.F.K. and J.F.R. developed the two-dimensional thermal model with magmatism. J.F.R. contributed a single-phase numerical code to compute the thermal structure of a subduction zone, to which D.R.J. added two-phase flow. R.F.K. and D.R.J. compared the model with petrological and heat flow observations. D.R.J. wrote the manuscript with R.F.K., and discussed the manuscript with M.T. and J.F.R. All authors jointly discussed and analyzed the data, results, conclusions, and implications.

Acknowledgements

The authors thank D. McKenzie, P. England, B. Hacker and J. Ague for comments on an earlier version of this manuscript. We would like to thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme Melt in the Mantle which was supported by EPSRC Grant Number EP/K032208/1. D.R.J. acknowledges research funding through the NERC Consortium grant NE/M000427/1 and NERC Standard grant NE/I026995/1. The research of R.F.K. leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement number 279925. M.T. received research funding from the Royal Society Newton International Fellowship. J.F.R. thanks the Leverhulme Trust for support. The authors would also like to thank the Deep Carbon Observatory of the Sloan Foundation.

S1 Petrological model of hydrous flux melting

In the 1D column model, we use a simple petrological model of hydrous flux melting, which is the dominant form of melting in subduction zones. The model was inspired by previous studies Hirschmann1999III ; Katz2003 ; Keller2016 , and is developed as follows. First, we restrict attention to a ternary system. The three components should not be thought of as identifiable minerals or oxides but rather as idealized components chosen to capture the physics in which we are interested. We start with two components that can be considered ‘refractory’ and ‘fertile’ Ribe1985 ; Hewitt2010 . To this system, we add a third component to represent volatiles. We initially take this component to be ‘water’ and we consider that the concentration of ‘water’ is relatively small. One role of this third hydrous component is to depress the solidus temperature.

Our second simplification is to use a quasi-linear phase diagram. This can be thought of as a linearization of the ternary phase loops used by Ref. Keller2016 about some initial composition.

Our third simplification is that the melting/solidification reactions happen sufficiently rapidly that a partially molten region is at local thermodynamic equilibrium. This implies that compositions of the coexisting solid and liquid phases are given exactly by the phase diagram.

S1.1 Mathematical description of phase diagram

The solidus temperature increases with increasing pressure at a rate $\gamma$. We linearize the dependence of the solidus on chemical composition. Since the sum of the concentrations of the components is unity, we need only specify two linear coefficients $M_{2},M_{3}$ for the fertile and water components respectively, both of which lower the solidus temperature. Thus the solidus temperature

[TABLE]

This expression can be rearranged to give, for example, the solidus concentration $c^{s}_{2}$ as a function of temperature, depth, and concentration of the third component. An interpretation of equation (S5) can be made by identifying $T_{s0}-\rho gz/\gamma$ with the solidus temperature of the refractory component at given depth $z$, which in this section is negative.

We assume that the liquidus concentration is related to the solidus concentration as follows:

[TABLE]

For the simplest case we take $\Delta c_{2,3}$ to be constants, but we will also consider generalizations.

S1.2 Choice of parameter values and implications for melting

We choose parameters in our model to constrain the degree of so-called ‘batch melting’ as a function of temperature and pressure:

[TABLE]

Batch melting refers to the degree of melting experienced by a sample raised to given temperature and pressure conditions assuming no extraction of melt. The composition $\widetilde{c_{0j}}$ is the composition of the solid mantle before the onset of melting. We then combine equations (S5)–(S7), which apply for each $j$, with equation (S8) to obtain

[TABLE]

A key quantity is the isobaric productivity $\partial F/\partial T$. If $\Delta c_{2}$ and $\Delta c_{3}$ are constants, then the isobaric productivity is a constant

[TABLE]

Thus melt is produced at a constant rate with increasing temperature. Linear models of two component melting already include this effect (e.g. Hewitt2010, ). It is worth noting that the productivity is reduced by the third, hydrous component.

In this formulation, volatiles do indeed depress the solidus temperature. However, in addition to depressing the solidus, volatiles are also associated with a ‘low-productivity tail’ Hirschmann1999III . The initial melting above the solidus temperature is less productive than later melting:

[TABLE]

The purely linear model does not satisfy this constraint, because the productivity is constant. Therefore, we generalize our model to allow for a low-productivity tail. Volatiles are incompatible and partition into the melt with a partitioning coefficient $D$ defined by

[TABLE]

where $D\ll 1$ for incompatible, volatile elements. We assume that $D$ is constant. However, $\Delta c_{3}$ is no longer constant, instead depending on composition, and hence pressure and temperature. Upon a little rearrangement, we find

[TABLE]

which can be rearranged to give a quadratic equation for $c_{3}^{s}$, recalling that $\Delta c_{3}$ is proportional to $c_{3}^{s}$. The degree of melting $F$ is no longer a linear function (however it can be computed explicitly using the quadratic formula so there is no computational difficulty, unlike more complex nonlinearities where iterative methods are required to solve for $F$). We can calculate the isobaric productivity at $F=0$ and $F=1$ and find

[TABLE]

We can interpret the effective isobaric productivity of the mixture as the harmonic mean of productivities associated with the fertile and volatile component. Typically, the contribution of the volatile component dominates at small $F$ and the fertile component dominates at large $F$ because

[TABLE]

Our formulation thus achieves the low-productivity tail expected physically; and it recovers the fertile–refractory system in the absence of volatiles.

To summarize, the degree of melting increases over the temperature range

[TABLE]

The gradient of the function $F(T)$ at these endpoints is given by equations (S14) and (S15).

Thus a limited number of parameters can describe a significant range of realistic melting behaviours, as shown in Fig. S1. For the anhydrous, fertile–refractory part of the system, we use $T_{s0}=1100$ ∘C, $\rho g/\gamma=4.5\times 10^{-3}$ ∘C/m, $M_{2}=700^{\circ}$C, $\Delta c_{2}=0.6$, and $\widetilde{c_{02}}=0.15$. This ensures that we match the anhydrous melting curve of Ref. Katz2003 , particularly around 3 GPa. For the hydrous part of the system, we use $M_{3}=2\times 10^{5}$ ∘C, $D=0.01$. These parameter values were chosen to roughly match the hydrous melting curve of Ref. Katz2003 , particularly around 3 GPa with 0.5 wt% water. Although the precise parameter values are in the right region for consistency with previous studies and their associated experimental libraries, as well as inferences from field observations, this parameterization is too simple to reproduce all the features observed experimentally. However, it can reproduce the two main features: solidus depression and a low-productivity tail.

S1.3 Generalized model: accounting for saturation in water

The addition of more water does not indefinitely lower the solidus, because eventually water becomes saturated in the liquid phase. The amount of water that dissolves increases with pressure; Ref. Katz2003 uses the formula

[TABLE]

where the pressure $P$ is measured in GPa. This is well constrained by experiment below 2 GPa, and constrained indirectly at higher pressures. This corresponds to a critical degree of melting and critical temperature below which the degree of melting drops rapidly to zero, as shown in Fig. S2c, for example.

Our modelling approach is to mimic this behaviour by modifying the phase diagram. We first compute the corresponding critical solid saturation point $c^{s}_{\mathrm{sat}}$, using equation (S19) for the liquid saturation and the partition coefficient of equation (S12). For temperatures below this point, we change the freezing point depression coefficient:

[TABLE]

where $M_{4}\leq M_{3}$. Note that the previous model is a special case $M_{4}=M_{3}$, and a eutectic-like phase diagram can be obtained by the special case $M_{4}=0$. In practice, we find $M_{4}=M_{3}/50$ makes a decent approximation to Ref. Katz2003 , as shown in Fig. S2. This means that the initial productivity near $F=0$ is a factor $M_{3}/M_{4}=50$ times greater. We use this generalized model to assess the significance of the increased productivity near water saturation in Fig. 3 of the main text.

S2 Further details of one-dimensional column model

In the context of a one-dimensional melting model, mass conservation imposes a strong constraint on the model behaviour in steady state. We adopt an extended Boussinesq approximation in which density differences between the phases are neglected except for their role in driving buoyant liquid segregation. There are several equivalent ways to present the following equations; we approach the problem by considering conservation in the liquid phase and in the two-phase composite.

Mass conservation gives

[TABLE]

where $\overline{x}=x^{s}(1-\phi)+x^{l}\phi$ denotes an average over the solid and liquid phases, with volume fractions $(1-\phi)$ and $\phi$ respectively. The vertical velocity is $w$, volumetric melting rate is $\Gamma$ and density is $\rho$. We first integrate equation (S22) to obtain

[TABLE]

where ${\overline{W_{0}}}$ is the total volume flux at the bottom of the melting column (which is not the motion of the solid phase alone, unlike in upwelling mantle columns used in the context of mid-ocean ridge magmatism). We follow the approach of Ref. Ribe1985 and define the quantity $F={\phi w^{l}}/{\overline{W_{0}}}$. Thus the scaled, liquid-phase volume flux is $F$ and the scaled solid phase volume flux is $(1-F)$.

We can recover our previous definition of $F$ in equation (S8) by considering conservation of species mass. For each component $j=1,2,3$,

[TABLE]

Note that, by summing equation (S24) over $j$ and comparing with equation (S21), $\sum_{j}\Gamma_{j}=\Gamma$. We integrate equation (S25) and use equation (S23) to obtain

[TABLE]

We then determine the degree of melting $F$, which is controlled by an energy equation and our phase diagram. One unusual feature of subduction zones is the non-monotonic temperature profile, which is largely controlled by the flow of the solid mantle. As described in the main text, we use a steady energy balance for a one-dimensional column

[TABLE]

with a volumetric source term $\rho c_{p}\Psi$ that represents the heat supplied by large-scale mantle corner flow. In the absence of melting, the final pair of terms on the right-hand side establishes a non-monotonic temperature profile. Note that we do not solve momentum equations because only the two-phase average velocity $\overline{W_{0}}$, which is constant as a result of mass conservation, appears in the heat equation (S27).

Next we observe that $\overline{W_{0}}F^{\prime}=\Gamma/\rho$ and $F^{\prime}$ is proportional to the isobaric productivity discussed previously, namely $F^{\prime}=(T^{\prime}+\rho g/\gamma)\partial F/\partial T$. We can better understand the system by rescaling the energy equation. We scale lengths by $H$ (the depth of the melting column), and the source term by $\kappa/H^{2}$. The dimensionless parameters involved are a Péclet number $\textrm{Pe}\,=H\overline{W_{0}}/\kappa$, a Stefan number $\textrm{St}\,=(L/c_{p})\partial F/\partial T$, a temperature change $\Delta T_{H}=\rho gH/\gamma$. Then the energy equation is

[TABLE]

A scaled version of the melting rate is

[TABLE]

which has units of degrees Kelvin. Equation (S28) is subject to boundary conditions on $T$ at $z=0$ and $z=1$. In general, the Péclet number is fixed but the Stefan number depends on temperature and pressure (hence depth), as well as the compositional parameters of our melting model. We plot the Stefan number in Fig. S3.

Our column-model approach is as follows. Extract a vertical temperature profile $T_{\mathrm{ref.}}$ from a single-phase mantle flow and thermal model of a subduction zone, as shown in the main article. We then calculate $\Psi=-T^{\prime\prime}_{\mathrm{ref.}}(z)$. To investigate the effect of melting, we solve the rescaled energy equation, focussing on the effect of varying the Péclet number and Stefan number (since the decompression term $\Delta T_{H}$ is well known). We present results in the main article.

S3 Further details of two-dimensional thermal model

Sources of fluids in subduction zones that trigger silicic magmatism are believed to be localized to particular depth ranges, associated with particular dehydration reactions in the subducting slab. Thus, in addition to the calculations presented in the main article, we also take three Gaussian magma flow profiles above the locations of the major dehydration reactions of the slab, with a position, magnitude and width suggested by Ref. Wilson2014 . We also consider the effect of all three sources combined.

As in the calculations in the main text, the principal result is that advective transport by magma substantially alters the thermal structure of subduction zones, as shown in Fig. S4. Flow associated with the peridotite source (a) is the most thermally significant, raising temperatures by over 200 K. Flow associated with the MORB source (b) raises temperatures near the trench by about 40 K; the gabbro source (c) is thermally insignificant. The peridotite source is strongest because it is associated with the largest magma flux. The MORB source is weaker because the flux is smaller and also because it occurs nearer the trench than the peridotite source, which means that the mantle wedge above the MORB source is slightly cooler. The gabbro source is especially weak because the flux is smaller, and because it is narrower than the other sources, and so tends to diffuse laterally more strongly. The combined set of sources (d) is dominated by the peridotite source, although there are also slightly elevated temperatures in the fore-arc region associated with the MORB source.