Articles | Volume 13, issue 9
Solid Earth, 13, 1415–1430, 2022
Solid Earth, 13, 1415–1430, 2022
Research article
05 Sep 2022
Research article | 05 Sep 2022

Assessing the role of thermal disequilibrium in the evolution of the lithosphere–asthenosphere boundary: an idealized model of heat exchange during channelized melt transport

Assessing the role of thermal disequilibrium in the evolution of the lithosphere–asthenosphere boundary: an idealized model of heat exchange during channelized melt transport
Mousumi Roy Mousumi Roy
  • Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87106, USA

Correspondence: Mousumi Roy (


This study explores how the continental lithospheric mantle (CLM) may be heated during channelized melt transport when there is thermal disequilibrium between (melt-rich) channels and surrounding (melt-poor) regions. Specifically, I explore the role of disequilibrium heat exchange in weakening and destabilizing the lithosphere from beneath as melts infiltrate into the lithosphere–asthenosphere boundary (LAB) in intraplate continental settings. During equilibration, hotter-than-ambient melts would be expected to heat the surrounding CLM, but we lack an understanding of the expected spatiotemporal scales and how these depend on channel geometries, infiltration duration, and transport rates. This study assesses the role of heat exchange between migrating material in melt-rich channels and their surroundings in the limit where advective effects are larger than diffusive heat transfer (Péclet numbers > 10). I utilize a 1D advection–diffusion model that includes thermal exchange between melt-rich channels and the surrounding melt-poor region, parameterized by the volume fraction of channels (ϕ), average relative velocity (vchannel) between material inside and outside of channels, channel spacing (d), and timescale of episodic or repeated melt infiltration (τ). The results suggest the following: (1) during episodic infiltration of hotter-than-ambient melt, a steady-state thermal reworking zone (TRZ) associated with spatiotemporally varying disequilibrium heat exchange forms at the LAB. (2) The TRZ grows by the transient migration of a disequilibrium-heating front at a material-dependent velocity, reaching a maximum steady-state width δ proportional to ϕvchannel(τ/d)n, where n≈2 for periodic thermal perturbations and n≈1 for a single finite-duration thermal pulse. For geologically reasonable model parameters, the spatiotemporal scales associated with establishment of the TRZ are comparable with those inferred for the migration of the LAB based on geologic observations within continental intra-plate settings, such as the western US. The results of this study suggest that, for channelized transport speeds of vchannel=1 m yr−1, channel spacings d≈102 m, and timescales of episodic melt infiltration τ≈101 kyr, the steady-state width of the TRZ in the lowermost CLM is ≈10 km. (3) Within the TRZ, disequilibrium heat exchange may contribute 10-5 W m−3 to the LAB heat budget.

1 Scientific motivation

During its long residence at Earth's surface, continental lithosphere is shaped by tectonic events such as rifting (including supercontinent break-up) and plate collision, undergoing profound changes in its physical and chemical state. In some cases, previously stable (undeforming) portions of the continental lithosphere may be destabilized. Based on the close association of magma infiltration with these events there is growing speculation that, under certain circumstances, melt–rock interaction may somehow weaken and perturb the stability of continental lithosphere (Hopper et al.2020; Wenker and Beaumont2017; Plank and Forsyth2016; Roy et al.2016; Wang et al.2015; Menzies et al.2007; Carlson et al.2004; Gao et al.2002; O'Reilly et al.2001). Recent work in subduction settings suggests that heat advection by magma transport into the overriding lithosphere is a fundamental process that determines the thermal structure of arcs and possibly aids thinning of the overriding plate (England and Katz2010; Perrin et al.2016; ReesJones et al.2018).

Fundamentally, this melt-enhanced weakening of continents is intertwined with the notion of thermal and chemical disequilibrium between infiltrating melts and the surrounding material, and therefore with transient processes that drive the lithosphere from one stable state to another. These processes, however, remain elusive. This paper explores one aspect of the problem, namely, thermal disequilibrium during infiltration of hotter-than-ambient melts into the base of the lithosphere, as a means of weakening and shaping the continental lithosphere from beneath. Specifically, I explore thermal disequilibrium between melt-rich channels and surrounding (melt-poor) material as a process to heat and modify the continental lithospheric mantle (CLM). The models explored here are a useful way to place constraints on the melt-transport scenarios for which a significant degree of disequilibrium heat exchange (driven by temperature differences between materials within and surrounding channels) may be important in the CLM at or near the lithosphere–asthenosphere boundary (LAB).

This study is inspired by evidence for the role of thermal disequilibrium from detailed field-based, petrologic, and geochemical studies in the Lherz and Ronda peridotite massifs (Bodinier et al.2008, 1990; LeRoux et al.2007; Lenoir et al.2001; Saal et al.2001; Vauchez and Garrido2001; LeRoux et al.2008, 2007; Soustelle et al.2009). Important conclusions from studies in the Lherz and Ronda peridotite massifs include the following: (1) “lherzolite” (named after its type-section in the Lherz massif), commonly regarded as pristine, fertile sub-continental lithospheric mantle, may be derived from refertilization of a depleted, harzburigitic parent (e.g., LeRoux et al.2007, 2008); (2) careful microstructural, geochemical, and petrologic work has documented the dominant effect of a steep thermal gradient associated with the region of contact and interaction between partial-melt-rich regions and the surrounding lithosphere. These workers provide a quantitative estimate of a transient thermal gradient (230C km−1, or more than an order of magnitude larger than a typical equilibrium geothermal gradient expected at the LAB). Indeed, the authors recognize this as a transient LAB and coin the term “asthenospherization” for disequilibrium heat exchange processes modifying the LAB. The spatial scale over which this disequilibrium heating is observed in Ronda (∼1 km) guides the mesoscale modeling approach here.

Additionally, this work is motivated by observations from the western US, which has undergone extensive magma-infiltration in Cenozoic time. Pressures and temperatures of last equilibration of Cenozoic basalts consistently point to depths that are at or below the present LAB (Plank and Forsyth2016), suggesting that melt transport from those depths upward through the lower CLM occurs in thermal disequilibrium. In the Big Pine volcanic field, for example, the inferred depth of the LAB decreases by >10 km in a timespan of <1 Myr (Plank and Forsyth2016), suggesting that the processes associated with this migration may be transient (LAB vertical migration rates in excess of ≈10 km Myr−1). More recently, Cenozoic melt- or fluid-enhanced thinning of the CLM in the western US has also been inferred from geochemical and isotopic data from volcanic rocks (Farmer et al.2020). Such processes of thermally driven erosion and migration of the LAB are also inferred in arc settings (e.g., Perrin et al.2016; ReesJones et al.2018). Motivated by these observations, a primary goal of this work is to quantify the role of transient, disequilibrium heating by infiltrating channelized melt as a mechanism for modifying the LAB and the lowermost CLM.

It has been argued that a permeability contrast (e.g., Holtzman and Kendall2010) or a change in magma mobility (e.g., Sakamaki et al.2013) across the LAB is likely to drive melts to pond and possibly drive the upward propagation of dikes that may freeze and heat the CLM (e.g., Havlin et al.2013; ReesJones et al.2018). Not all infiltrating melts would freeze, however, and some component of hotter-than-ambient melts may be transported in thermal and chemical disequilibrium into the CLM via established channels or pathways (e.g., LeRoux et al.2007; Schmeling et al.2018). Thermal disequilibrium during melt transport is expected to become important within the CLM as the degree of channelization and the relative melt–solid velocity increases (e.g., Schmeling et al.2018; Chevalier and Schmeling2022). In this work, I am not concerned with the emergence and development of these channel networks within the lower CLM, nor the deeper processes that transport melt from a sub-lithospheric melt-generation zone (e.g., Aharonov et al.1995). Instead, the starting point of this study is the observation that high-porosity, melt-rich channels are an important part of melt–rock interaction in the CLM (e.g., Soustelle et al.2009; LeRoux et al.2007). An important limitation of the models, therefore, is that channelization is imposed via parameters that control a heat transfer coefficient. The simplicity of the models, however, allows us to focus on the implications of significant thermal gradients between melt-rich channels and their surroundings. Although others have also argued for the important role of thermal disequilibrium in melt–rock interaction (Chevalier and Schmeling2022; Keller and Suckale2019; Wallner and Schmeling2016; Schmeling et al.2018), this study provides a quantification of the role of thermal disequilibrium at the LAB based on observational constraints discussed above. This work builds on the 1D model in Roy (2020) (which did not consider axial conduction) and includes both a thermal contrast that drives heat exchange and axial diffusion terms (e.g., following Chevalier and Schmeling2022) in order to explore thermal equilibration over long timescales of melt infiltration (>103-4 years) into the lowermost 1–10 km of the CLM (e.g., stage 3, large Péclet numbers in Chevalier and Schmeling2022).

The 1D model abstracts the complex geometry of the melt–rock interface and therefore differs from other descriptions of disequilibrium heat exchange (Wallner and Schmeling2016; Schmeling et al.2018; Keller and Suckale2019). Similar to chemical transport models that use a linear driving term for chemical disequilibrium (e.g., Hauri1997; Kenyon1990; Bo et al.2018), the models below assume a linear thermal driving term (e.g., Schumann1929; Kuznetsov1994; Spiga and Spiga1981). In this study, however, only the role of thermal disequilibrium is included and I ignore important processes such as chemical disequilibrium (e.g., Hauri1997; Kenyon1990; Bo et al.2018). The basic results of the 1D model are presented below, followed by further discussion of their limitations and implications. Although idealized, the models provide a first-order assessment of the temporal and spatial scales over which thermal disequilibrium can play a role in warming and therefore weakening the lowermost CLM. A key finding is that the lowermost portion of the CLM may evolve into a long-term thermal reworking zone (TRZ) driven by disequilibrium heat exchange in channelized melt transport. The factors that determine the rate at which the TRZ is formed and its spatial scale are estimated from the models and are compared to geologic observations within the western US, specifically geochemical and petrologic evidence for the upward migration of the LAB during Cenozoic melt–rock interaction (Plank and Forsyth2016; Farmer et al.2020).

2 Model of disequilibrium heat transport

The starting point for the models explored here is a simple, 1D theory of heat exchange in packed porous beds by Schumann (1929), building on Roy (2020), where fluid moves within the pores of a matrix of solid grains, and the only variability in temperature (and temperature contrast between solid and fluid) is in the transport direction. In Schumann (1929), the thermal evolution of the system is governed mainly by heat exchange across the solid–fluid interfacial surface. In this study, I present a re-interpretation of Schumann's model and of the effective heat transfer coefficient. Instead of considering fluid moving in pores between solid grains, the system of equations from Schumann (1929) may be used to describe thermal disequilibrium between material within high-porosity channels and outside channels. In other words, here “fluid” is interpreted to be in-channel material and “solid” is material outside channels (for simplicity, however, I retain the subscripts f (in-channel) and s (outside channels)). The goal here is to describe the relative importance of advective transport over length scales that are comparable to the channel spacings, so we define a transport velocity, vchannel, as an average relative velocity across channel walls. As described below, this reinterpretation also extends to the physical meaning of the effective heat transfer coefficient, where now the geometry across which the transfer occurs must take into account the channel geometry and spacing. Furthermore, following Schumann (1929), heat exchange is assumed to be linearly proportional to the local temperature difference between solid and fluid (see also Roy2020). Unlike Schumann (1929) and Roy (2020), however, this study includes thermal diffusive effects to account for (axial) conduction within channels and within the material outside of the channels. These arguments lead to coupled equations for the temperature outside channels, Ts, and within channels, Tf:


where vchannel is the average transport velocity of material within melt-rich channels relative to the melt-poor surrounding material, ϕ is a fluid volume fraction, λf and λs are thermal conductivities, cf and cs are the heat capacities per unit volume (heat capacitances) at constant pressure (cf=cpfρf and cs=cpsρs), and x is the position coordinate in the transport direction. Note that the geometry of the solid–fluid interface is not treated in detail but is idealized in the channel volume fraction, ϕ, the channel spacing, d, both of which control an effective (across-channel-wall) heat transfer coefficient, K, discussed below.

One advantage of “coarse-graining” in the Schumann (1929) model (from the pore scale to macroscopic channels) is its simplicity and that it has been investigated in numerous previous studies. There exist analytic solutions for Eqs. (1) and (2) in limiting cases, particularly for large Péclet numbers with axial diffusion terms ignored (without the last terms in Eqs. 1 and 2) (Spiga and Spiga1981; Kuznetsov1994, 1995a, b, 1996). This re-interpretation of the model must also be accompanied by an appropriate reinterpretation of the heat transfer coefficient, K, made possible because in the framework above the geometry of the interfacial surface is not explicitly specified. The reinterpreted model is applied to a semi-infinite domain where fluid transport occurs in high-porosity channels aligned in one direction (Fig. 1). The high-porosity channels are assumed to occupy a constant volume fraction, ϕ, within which material moves with a constant (average) velocity vchannel relative to the surrounding stationary material outside the channel (volume fraction 1−ϕ). The model domain may be thought of as co-moving with the reference frame of material outside the channels. Because of the assumptions built into the 1D approach, the results below are applicable to physical situations where transport is dominantly in the along-channel direction and any motion of material outside channels is steady. The model assumes that the average channel geometry is unchanging within the domain and that the material outside the channels is initially at a uniform temperature. The models and their interpretations, therefore, are used here to assess the role of thermal disequilibrium as melts are transported (≈10–20 km or so) within the lowermost CLM above the LAB within intraplate settings.

Figure 1Illustration of the 1D model specified by (1) material parameters within and outside of channels (heat capacities cpf and cps, thermal conductivities λf and λs, and densities ρf and ρs), (2) average in-channel velocity vchannel, and (3) geometric parameters such as channel volume fraction ϕ and channel spacing d. The effective heat transfer coefficient K is a function of heat capacitances (cpf/ρf, cps/ρs) ϕ and d, scaling as d−2 (see Sect. 2.1); large K corresponds to large channel wall area per unit volume (e.g., small d) and vice versa. The input in-channel temperature vs. time functions considered in this study are shown in the three graphs below: (I) step function, (II) sinusoid, (III) finite-duration pulse.


In this model, the terms involving the thermal contrast (TfTs) in Eqs. (1) and (2) represent heat transfer across the walls of the channels, between material within and outside the channels. This heat exchange depends on material parameters that govern thermal diffusion perpendicular to the transport direction, and on the geometry of the channels themselves (e.g., sinuosity, spacing). The detailed geometry of the channel walls (the relevant interfacial surface here) is not specified in this simple 1D treatment, but is parameterized by the heat transfer coefficient, K (Fig. 1). Therefore, K is a proxy for the efficiency of conduction perpendicular to the transport direction, and the geometry of the channel wall interface, namely the wall area per unit volume, controlled by the spatial scale of channelization, d (see Sect. 2.1).

As illustrated in Fig. 1, a large value of K may represent efficient heat exchange as in the case of many channels separated by a small distance. Conversely, a low value of K would represent inefficient exchange, as in the case of a larger characteristic separation between the channels. In the following sections, I consider the physical meaning of K and also present a non-dimensionalization of the Eqs. (1) and (2) based on characteristic length and timescales in the problem. The coefficients of the temperature-contrast terms on the right hand sides of Eqs. (1) and (2) specify the timescales of heat exchange within channels, tf=1/Kf=ϕcf/K, and outside channels, ts=1/Ks=(1-ϕ)cs/K. Additionally, a characteristic length scale emerges out of the relative motion across channel walls, vchannel/Kf. These characteristic length and timescales are used to non-dimensionalize Eqs. (1) and (2) (see Sect. 2.2) and obtain the results presented below.

2.1 Heat transfer coefficient

In this section I consider the meaning of the heat transfer coefficient, K, and the related constants, Ks=K/(cs(1-ϕ)) and Kf=K/(cfϕ) in Eqs. (1) and (2). Note that Kf and Ks have dimensions of inverse time, and they both depend on the heat transfer coefficient K. Physically, K represents the amount of heat transferred across channel walls per unit time, per unit volume, and per unit difference in temperature (in  Schumann1929, this exchange is across the solid–fluid interface). The factors that determine K can be illustrated by considering that the heat transfer rate across channel walls must depend on the geometry of walls and also on the effective thermal conductivity of the channelized domain.

Although the geometry of the channels may be complex, I consider one aspect of it, namely, the specific wall surface area (wall area per unit volume), asf, which is a function of the spatial scale of channelization. In the grain-scale porous flow case considered in Schumann (1929) for example, if the solid matrix is made of spheres with an average particle diameter p, then the specific area for a grain is S0=6/p, so asf=S0(1-ϕ)=6(1-ϕ)/p (Dullien1979). This sets a limit for channels with channel spacing d, where we shall assume that the specific surface area is asfA(1-ϕ)/d, where A is a number that is between 2 (for planar sheet-like channels with small volume fraction ϕ) and 6 for transport around spherical regions (Dixon and Cresswell1979; Schmeling et al.2018). Whereas the specific wall area is a geometric factor, the effective conductivity of the medium depends on the Nusselt number, Nu. Theoretical arguments in Dixon and Cresswell (1979) show that the effective thermal conductivity may be written in terms of the individual in-channel and out-of-channel thermal conductivities λf and λs (equivalent to considering the channels and non-channel regions in parallel):

(3) 1 C eff = 1 N u λ f + 1 β λ s ,

where β=10 for spherical matrix grains, 8 for cylinders, and 6 for slabs (Dixon and Cresswell1979). Therefore, the range of β=6 to 10 represents the highly channelized vs. porous flow end-members. For slow flows (Reynolds number Re≪100), Handley and Heggs (1968) argue that Nu ranges from 0.1 to 12.4 (Dixon and Cresswell1979) (Table 1). The relevant quantity that determines K is an effective “conductance” Ceff/d, so that K=Ceffasf/d,

(4) K = 1 N u λ f + 1 β λ s - 1 A ( 1 - ϕ ) d 2 ,

a product of a material-dependent quantity and a geometry-dependent quantity.

Lesher and Spera (2015)Lesher and Spera (2015)Lesher and Spera (2015)Lesher and Spera (2015)Lesher and Spera (2015)Pec et al. (2017)Rutherford (2008)(Handley and Heggs1968)Dixon and Cresswell (1979)Dullien (1979)LeRoux et al. (2008)

Table 1Material properties and constants used in calculations

Download Print Version | Download XLSX

Turning now to physical properties relevant to the transport of melts through the lithosphere, the channelized domain may be thought of as consisting of a mixture of melt+grains throughout, but with variable grain-scale porosity φ. Specifically, the fluid-rich channels would have a higher φf than the surroundings φs. The thermal conductivity inside and outside the channels would then be a volume average in each, e.g., inside channels, λf=φfλmelt+(1-φf)λgrain, where λmelt and λgrain are the values for, say, basaltic melt and peridotitic grains. Similarly, outside channels the volume average would be λs=φsλmelt+(1-φs)λgrain. Here I do not specify reasonable ranges of φf and φs (in general, they will both be within [0,1]), but rather I focus on determining an upper limit to the role of disequilibrium heat exchange across channel walls. Therefore, to explore the end-member case, I take φf=1 and φs=0, so that λf=λmelt and λs=λgrain. Using a reasonable conductivity for basaltic magma of λmelt=1 W m−1 K−1 (Lesher and Spera2015), λgrain=2.5 W m−1 K−1) for the solid grains, and taking Nu=0.1 to 12.4, β=6, 8, or 10, and A=6, we find that the effective conductivity Ceff, asf, and K are within the ranges shown in Fig. 2. (Note that choosing A≈2 would change K by less than an order of magnitude.) As suggested by Eq. (4), Kd-2 and strongly decreases with increasing spatial scale of channels (Fig. 2c). Although Fig. 2 theoretically explores the full range of parameters and their effect on K, in practice these parameters are not independent and should be chosen based on the geometry considered (e.g., Chevalier and Schmeling2022). In the models in Sect. 3 I set β=6 and A=2 as suggested for 1D channels (Dixon and Cresswell1979; Schmeling et al.2018). As discussed in Schmeling et al. (2018), the heat transfer coefficient should depend not only on d the channel spacing, but also a length scale set by the thickness of a microscopic thermal boundary layer at fluid–solid interfaces, specifically, KCeff/ (a boundary layer dimension) (Schmeling et al.2018). This boundary layer thickness is a function of time and only at timescales that are long relative to a characteristic thermal response time will the boundary layer encompass the entire region between channels Schmeling et al. (2018). Therefore, by taking Kd-2, our models must be confined to timescales that are large relative to the material-dependent thermal response timescale (e.g., 1/Ks), a requirement that is met in all thermal perturbations considered in this study.

Figure 2(a) Effective thermal conductivity, Ceff in Eq. (3), as a function of Nusselt number; (b) geometric factor, asf, as a function of channelization scale, d; and (c) heat transfer coefficient, K, as a function of channelization scale d. For a fixed d, the dashed lines in (c) delineate the variation in K for the range of β values in (a) and ϕ values in (b), illustrating that K is mainly controlled by d, rather than the other parameters. Gray shading indicates the range of channel spacings relevant to models with Pe>10 presented in this study.


To decide on a range of K values appropriate to the LAB, I turn to geologic observations of the scale of channelization in exhumed portions of the lower CLM. Structural, petrologic, and geochemical data from the Lherz Massif suggest that melt–rock interaction has driven refertilization of a harzburgite body into lherzolite (LeRoux et al.2007, 2008). In the field, the lherzolite bodies are separated from each other by distances of several tens of meters and this is also the spatial scale of isotopic disequilibrium between metasomatizing fluids and the harzburgite parent material (LeRoux et al.2008). With this as a proxy for the spatial separation of fluid-rich channels, I choose a broad range for the relevant spatial scale of channelization, d=100 to 103 m (1 m to 1 km channel spacing). The corresponding range of the heat transfer coefficient in the models is therefore K10-5 (large d) to 101 (low d) in W m−3 K−1 (Fig. 2). In the following, material properties, channel volume fraction, channel velocity, and heat transfer coefficient are fixed for each calculation (Table 1), but we confine d>102 m in order to ensure an effective Péclet number Pe>10. Taking typical parameters ϕ=0.1, d=1000 m, and corresponding K values in Fig. 2, typical timescales of response are ts=1/Ks9800 years and tf=1/Kf1200 years, short compared to the timescales of geologic events and the timescale of sinusoidal and pulse-like thermal perturbations considered here.

2.2 Non-dimensional system

To non-dimensionalize the system of Eqs. (1) and (2) (following Spiga and Spiga1981), we define the normalized relative temperature, Tf=(Tf-T0)/ΔT and Ts=(Ts-T0)/ΔT, where T0 is reference temperature and ΔT is a temperature perturbation (described below). We also introduce the dimensionless position, x=xKf/vchannel, a dimensionless time, t=Kst, and the weighted heat capacitance ratio, ζ=Ks/Kf=ϕcf/(1-ϕ)cs (from Eq. (5), we see that 1/ζ is a non-dimensional velocity). The non-dimensional versions of Eqs. (1) and (2) are now Eqs. (5) and (6):


Analytic solutions for this set of equations have been derived for a number of limiting cases, particularly for Df=Ds=0 (Spiga and Spiga1981; Kuznetsov1994, 1995a, b, 1996), and were used to test the numerical calculations in this study.

The terms in square brackets in Eqs. (5) and (6) represent dimensionless coefficients that govern the diffusion terms, Df and Ds. Using the definition of ζ, we can further simplify the coefficient of the diffusion term Ds and express it in terms of Df as

(7) D s = D f 1 - ϕ ϕ λ s λ f .

It is clear that for a given temperature difference (Tf-Ts) the behavior of Eqs. (5) and (6) is governed by Df and ζ (1/ζ is a dimensionless in-channel velocity), thermal conductivities λs and λf, and the channel volume fraction, ϕ. We can define an effective Péclet number for the problem as the product of the velocity vchannel, the characteristic length scale vchannel/Kf=vchannel(ϕcf)/K, divided by the thermal diffusivity of the channel material λf/cf,

(8) P e 1 = v channel 2 c f 2 ϕ λ f K = 1 / D f ,

where heat exchange due to thermal disequilibrium across channel walls (disequilibrium heating) will be important when Pe1≫1. Since K is a function of user-specified material and channel geometry parameters, we may further write Pe1 as

(9) P e 1 = v channel 2 c f 2 λ f 1 N u λ f + 1 β λ s d 2 ϕ A ( 1 - ϕ ) .

Alternatively, one may also define an effective Péclet number as the product of vchannel2, a characteristic timescale 1/Ks=cs(1-ϕ)/K, divided by the thermal diffusivity of material outside channels λs/cs,

(10) P e 2 = v channel 2 c s 2 λ s 1 N u λ f + 1 β λ s d 2 ϕ A .

These definitions are limited to consideration of channelized flow and more complex dependence between Pe, and ϕ is suggested in the case of tubes or pores (e.g., Chevalier and Schmeling2022). Using values of Nu=12, β=6, A=2, with ϕ and material properties in Table 1, for a given vchannel, d, Pe1 is ≈2 to 3 times larger than Pe2. Hereafter, I use the definition of Pe2 when referring to Péclet number Pe.

Note that the definitions in Eqs. (9) and (10) referred to as effective Péclet numbers, to distinguish them from the Péclet number one, can define for axial transport in the channels, Pechannel=vchanneldcf/λf. By using the characteristic times 1/Kf or 1/Ks, these definitions explicitly consider timescales associated with thermal exchange perpendicular to the transport direction, which depends on material parameters both inside and outside channels. In this manner, the definitions in Eqs. (8) to (10) take into account the key role played by heat exchange governed by K that works to reduce the thermal contrast across the channel walls.

2.3 Numerical method

Equations (5) and (6) are solved numerically in explicit time using a leap-frog method. The 1D domain x=[0,L] with L=102 to 103 is discretized typically with N=1000–5000 elements, and the maximum time for each run is tmax>Lζ (for thermal perturbations considered below that have a characteristic time τ, I use tmax>Lζ+2Ksτ). For a given element size dx=L/N, second-order finite differences are used for all spatial derivatives. The time step dt is chosen (empirically) to be small enough to avoid numerical dispersion (for Pe>1, it is sufficient to choose dt=dx×10-4). The code is written in MATLAB (R2021b), and results below are confirmed using different grid resolutions. Step-function perturbation models below without axial diffusion (Df=DS=0) were benchmarked (using leap-frog with upwind differencing for the spatial derivatives) against analytic solutions in Schumann (1929).

3 Results

In the models below, I consider a 1D domain (x≥0) that is initially in equilibrium Ts=Tf=0; material outside the channels is stationary, and only the in-channel material is moving, at speed vchannel. At x=0 and t=0+, the in-channel material is subjected to a thermal perturbation in Tf. The non-dimensional Eqs. (5), (6), and (7) are used to determine the transient thermal equilibration between material inside and outside channels, but we limit consideration to models where advective transport dominates and Pe>10. I consider transient thermal evolution with three scenarios of influx of hotter-than-ambient melt or fluid, in order of increasing complexity (Fig. 1): (I) a step-function increase in temperature of the in-channel material; (II) a sinusodal temperature perturbation; and (III) a finite-duration, constant-amplitude thermal pulse. The perturbation in Tf disturbs the initial steady-state condition in the domain starting at t=0+ as material with a perturbed temperature enters channels at x=0 with a positive perturbation amplitude ΔT (Fig. 1). The perturbation “front”, the farthest extent of channel material that entered at x=0 with perturbed temperature, is at xpert=vchannelt (or xpert=t/ζ; where the dimensionless velocity inside channels is 1/ζ; see Sect. 2.2). (Note that the location of maximum disequilibrium heat exchange (Tf-Ts)max) moves at a different speed, Vdiseqm, discussed below.) In each case, the physical, user-specified quantities that specify the model are the channel spacing d, channel volume fraction ϕ, the in- and out-channel thermal conductivities λf and λs, heat capacitances cf, and cs, and the speed of the in-channel material vchannel (see values in Table 1). Since the material properties cs and cf are held constant, there is a unique mapping between ζ and ϕ, the channel volume fraction (for ϕ=0.01 to 0.20, ζ=0.0096 to 0.2376; Table 1).

3.1 Response to a step function

The domain is initially at steady state in equilibrium at temperature Ts=Tf=T0 (Ts=Tf=0), and at t=0 the temperature of the in-channel material entering at the inflow x=0 is perturbed so that Tf(x=0,t>0)=T0+ΔT (or Tf=1). As material with higher temperature enters the channels, a transient thermal response occurs as the material inside and outside channels begin to equilibrate (material in channels cool, surroundings heat). A transient disequilibrium zone is observed to travel into the domain (Fig. 3a and b). Ahead of this disequilibrium zone, the channels are in equilibrium with the surroundings at the initial ambient temperature, Ts=Tf=0. Behind this zone, the channels are in equilibrium with the surroundings at the inlet temperature, Ts=Tf=1. The models compared in Fig. 3a and b (same ϕ but different channel spacing d) show that the transient disequilibrium zone is observed to travel into the domain at the same Vdiseqm, controlled by ϕ, but the observed broadening of initially steep thermal profiles during transport differs, suggesting it is controlled primarily by the channel spacing d, and therefore the effective heat transfer coefficient K. Although diffusion plays a role in the broadening of the thermal profiles downstream, it is important to note that, unlike a simple advection–diffusion equation (where at large Péclet number we might not expect as much broadening of an initially sharp pulse), both terms on the right hand sides of Eqs. (5) and (6) will drive broadening in this model. Even in the absence of diffusion, therefore, thermal contrast between the material inside and outside channels causes shallowing of steep thermal gradients during transport (e.g., Fig. 3c and d; see also Roy2020). At large Pe this disequilibrium heat exchange dominates over axial conduction.

Figure 3(a) Normalized temperature profiles, Ts (dashed lines) and Tf (solid lines), for models with in-channel velocity vchannel=1 m yr−1, at times t=50, 125, and 200 kyr after a step-function perturbation in Tf. The temperature profiles transition between the incoming channel material normalized temperature (=1, left) and the initial ambient temperature (=0, right). We compare temperature profiles for two cases with the same channel volume fraction ϕ=0.12 but different channel spacing, d, and therefore heat transfer coefficient K and Péclet number, indicated. The transition region (e.g., highlighted in gray for the d=500 m model at t=50 kyr), has width, w, that is larger for smaller K (large d) and increases over time. (b) Profile of the temperature difference Tf-Ts between channels and surroundings at times indicated, for the same calculations in (a). The degree of disequilibrium (red and blue arrows) is greater for smaller K (higher Pe) and decreases as a function of time. Panels (c) and (d) compare the results for d=500 m (red lines in (c) and (d) are identical to those in (a) and (b)), but using a modified model where the axial diffusion terms in Eqs. (5) and (6) are neglected (Df=Ds=0).


Following an initial lag time (when the maximum thermal contrast is at x=0), the disequilibrium zone (marked by the peak in the Tf-Ts function, (Tf-Ts)max in Fig. 4b) migrates inward migration at a steady speed Vdiseqm, a fixed fraction of vchannel that depends on ϕ, the channel volume fraction (Fig. 4a and b). The ratio Vdiseqm/vchannel varies linearly with ϕ (Fig. 4c). In the near-equilibrium limit (Tf-Ts0), Kuznetsov (1994) shows that the shape of the temperature difference function (Fig. 4b) approaches a Gaussian function with width that depends on t and the zone of disequilibrium migrates at speed vchannelcf/(ϕcf+(1-ϕ)cs). Our models show that, when there is significant disequilibrium, the zone of disequilibrium migrates with a rate given by Eq. (11) (Fig. 4c), independent of K, the heat transfer coefficient. Empirically, a key result is that the location(s) of maximum disequilibrium (Tf-Ts)max and therefore the greatest heat exchange progress inward into the domain at a rate given by material properties, the channel volume fraction, ϕ, and the in-channel velocity, vchannel,

(11) V diseqm v channel c f c s c f ϕ c f ϕ + ( 1 - ϕ ) c s .

Although Vdiseqm is independent of K, it is important to note that the degree of disequilibrium is not. Figure 3 illustrates the dependence on K for the specific case where the in-channel velocity vchannel=1 m yr−1, and d=500 to 1000 m, which correspond to K=4.7×10-5 and 1.2×10-5 W m−3 K−1, respectively. A second key result that emerges is that the degree of disequilibrium decreases exponentially as the zone of disequilibrium migrates inward. This spatial decay is observed in Fig. 3b; however, it is quantified below, considering periodic thermal perturbations.

Figure 4(a) Normalized temperature within channels (solid lines) and surroundings (dashed lines) at different times (indicated) following a step-function perturbation in Tf. The models shown here have the same in-channel velocity vchannel=1 m yr−1 and channel spacing d=500 and are compared at times t=50, 125, and 200 kyr as in Fig. 3. The two cases here differ, however, in the channel volume fraction, ϕ (indicated), illustrating that the migration rate of the zone of disequilibrium is a function of ϕ. (b) Normalized temperature difference Tf-Ts between channels and surroundings as a function of position, shown for the same times as in (a). (c) Normalized migration rate of the zone of disequilibrium as a function of channel volume fraction, ϕ. Red dot is for ϕ=0.12, corresponding to bright red lines shown in (a) and (b), and in Fig. 3.


3.2 Response to a sinusoidal thermal perturbation

Here we consider a second scenario where the fluid entering the domain is hotter than the ambient initial temperature, but the thermal contrast varies sinusoidally. Although this is an idealized condition, it may be interpreted to represent periodic pulses of high-temperature material entering into fluid- or melt-rich channels. Since any continuous time-varying thermal history at the inflow may be represented as a sum of sinusoids, this scenario also helps build intuition regarding the inherent length and timescales of equilibration. Sinusoidal thermal pulses introduce a new timescale into the problem, the period τ, and the relevant timescale to compare to is 1/Ks, the longest response timescale in the domain (on the order of 103 years; see Sect. 2.1), associated with the thermal response of the material outside channels.

Periodic thermal perturbations that might represent melt infiltration pulses lasting 104 to 106 years are characterized by a region of spatially varying temperatures: a thermal re-working zone (TRZ) (Fig. 5). A key result is that thermal pulses with periods that are long compared to 1/Ks penetrate farther into the domain than shorter period oscillations (Fig. 5). The non-dimensional period, τKs=τK/(1-ϕ)cs, controls the length scale, δ, over which thermal oscillations penetrate into the domain. The wavelength of these temperature oscillations is set by the period τ, λ=vchannelτ. The penetration distance of the oscillations, δ, is the maximum width of the TRZ. At short timescales, when tδ/Vdiseqm, there is one zone of disequilibrium (the TRZ), bounded by the inward-traveling edge of the region of spatially oscillating temperatures. The TRZ is initially narrow but widens at a rate Vdiseqm, eventually reaching a maximum width, δ (akin to a thermal “skin depth”), at time δ/Vdiseqm. When tδ/Vdiseqmt, there are two zones of disequilibrium: a stationary one bounded by the inlet with width δ (the TRZ, e.g., the wide region of gray at the left edge of the domain in Fig. 5a and b) and a migrating zone that continues inward at Vdiseqm (inset in Fig. 5b). The amplitude of the temperature oscillations decay with distance, but at each location in the TRZ the amplitude is constant, once oscillations are established (Fig. 5). This effect is identical to the exponentially decaying degree of chemical disequilibrium obtained using a very similar set of equations to Eqs. (2) and (1) in the analytic solutions of Kenyon (1990).

Figure 5(a) Normalized temperature profiles in-channel Tf (solid lines) and out-of-channel Ts (dashes), at times indicated, for a calculation with in-channel velocity vchannel=1 m yr−1, channel volume fraction ϕ=0.12, channel spacing d=500 m, and heat transfer coefficient K and Pe as indicated. For the chosen parameters, the response timescale is 1/Ks=ts2.6 kyr. Results are shown for two different sinusoidal thermal variations in the incoming in-channel material with (normalized) oscillation periods: τ/ts=20 (red, shown at t=350 kyr) and τ/ts=100 (blue, shown at t=700 kyr). The region of sinusoidal thermal profile (gray shading) has reached a steady-state width for the τ/ts=20 case but not for τ/ts=100. (b) Normalized temperature difference Tf-Ts between channels and surroundings as a function of position, shown for the same times as in (a). The thermal reworking zone (TRZ) (gray shaded region for τ/ts=20) has spatial oscillations in Tf-Ts, with amplitudes that decrease over a length scale δ, the width of the TRZ: δ is larger for longer period (blue) and shorter for shorter period (red). The degree of disequilibrium (Tf-Ts) is oscillatory in the TRZ, with decaying amplitude over width δ. Inset: Illustration of a steady-state TRZ (once it has widened to a width δ). During the widening phase, the TRZ width increases at a rate Vdiseqm. Once it has grown to maximum width, δ, beyond x=δ the domain is at equilibrium (Ts=Tf=0.5), except for a region that continues to migrate inward at Vdiseqm (arrows).


As we might expect, the maximum width of the TRZ, δ, is set controlled both by the non-dimensional oscillation period, τ/ts and by the heat transfer coefficient, K. In the absence of the axial diffusion terms, previous work (Spiga and Spiga1981; Kenyon1990) has shown that

(12) δ = c f ϕ v channel K ( τ / t s ) 2 4 π 2 ,

noting that Kd-2 (2.1) and ts=1/Ks=cs(1-ϕ)/K; the expression above suggests that, for fixed vchannel, δ is proportional to (τ/d)2 in the absence of axial conduction. With axial conduction, in the limit of large d (large heat transfer coefficient and Péclet number), this scaling is confirmed by the numerical results for sinusoidal periodic perturbations (square symbols in Fig. 7). A consequence of the inward-decaying degree of disequilibrium between material inside and outside of channels is that disequilibrium heat exchange occurs only within a certain distance, δ, of the inlet (x=0), within the TRZ (Fig. 5).

3.3 Finite pulse, duration τp, amplitude ΔT

Here we consider the fate of a finite-duration thermal pulse, representing episodic infiltration of melts that are hotter than the ambient CLM. This scenario introduces a timescale into the problem, the pulse duration τp, implemented here as a tanh-function (1/2)tanh(t-τp/2)w-tanh(t-3τp/2)w with a characteristic growth–decay timescale w that depends upon τp, e.g., w=0.1τp in these models. The results for a finite-duration pulse, exemplified in Fig. 6, are consistent with the characteristic transient behavior already apparent in the sine and step-function perturbations above. The thermal perturbation Tf in the channels at x=0 km (the inlet) is distorted as it proceeds into the domain, broadening in width and decaying in amplitude (Fig. 6a). Similarly, the zone of disequilibrium heat exchange (marked by (Tf-Ts)max) migrates into the domain at a rate given by Eq. (11) and decays during transport (Fig. 6b).

Figure 6Two views of the thermal response to a finite-duration perturbation (Tf(x=0,t) is a tanh-function) with duration τp=25 kyr, in a model with vchannel=1 m yr−1, channel volume fraction ϕ=0.12, and channel spacing d=500 m. (a) Normalized temperature outside channels (Ts, solid lines) and within channels (Tf, dashed lines) at various times as indicated. (b) The degree of disequilibrium, Tf-Ts as a function of position, the same times as in (a). Curved dashed line illustrates an exponentially decaying envelope that is used to estimate δ, the width of the TRZ. For clarity, intermediate temperature profiles used to estimate δ are shown in the inset. (c) Temperature–time history for a model with for a model with in-channel velocity. Normalized temperature outside channels (Ts, solid lines) and within channels (Tf, dashed lines) at varying locations within the domain (x, as indicated) are plotted as a function of time since contact with material that entered channels at x=0,t=0. (d) The degree of disequilibrium, Tf-Ts for the same x locations as in (c), plotted as a function of time since contact with material that entered channels at x=0,t=0.


As with the sinusoidal perturbation, these effects lead to a TRZ that widens over time to a maximum width, δ, that scales with the duration of the thermal pulse τ, channel spacing d, and volume fraction ϕ. Figure 7 illustrates how, for fixed vchannel=1 m yr−1, independent parameters d, ϕ, and τ lead to variable TRZ widths, δ. For a given channel spacing d, which strongly controls the heat transfer coefficient, Kd-2 (large d, small K), I consider a range of plausible channel volume fractions ϕ=0.01 to 0.2 (Table 1). Whereas δ is proportional to τ2 for large d in the pure sinusoid case (square symbols in Fig. 7), the finite-duration pulse results point to an n value likely less than 2, closer to 1 (Fig. 7). This is likely due to the fact that sinusoidal perturbations represent a higher and continuous energy input into the system compared to a finite-duration pulse, suggesting that in the case of multiple (episodic) pulses of melt infiltration, the scaling exponent is likely to be between 1 and 2. To obtain a TRZ width between 1 and 10 km, we require channel spacings around d=100 (with ϕ=0.1 to 0.2), d=500 m (with ϕ<0.16), or d=1000 m (with ϕ<0.02) and thermal pulse durations of 101 to 102 kyr (Fig. 7).

Figure 7Width of the TRZ, δ for calculations with vchannel=1 m yr−1 transport velocity and variable channel spacing, d. Numerically derived values of δ are obtained by fitting an exponential decay to the maximum degree of disequilibrium max(Tf-Ts), e.g., Fig. 6b. Results compiled here are for periodic thermal perturbations with period τp (squares) and for finite-duration thermal perturbations with duration τp (triangles). For clarity, models with fixed ϕ=0.12 and varying d are shown for the periodic perturbations. For finite-duration pulses, results shown include variable d and ϕ, as indicated. For given d and τp, δ is a function of channel volume fraction ϕ=0.01 to 0.2 (Table 1); the color saturation of the symbols corresponds to ϕ: darkest =ϕmin=0.01 and lightest =ϕmin=0.2. A few representative thin black dashed lines are drawn to connect models that have the same channel geometry parameters d and ϕ, but differ only by the incoming thermal pulse duration, τp. The slopes of the thin dashed lines therefore illustrate how δ scales with τp. Light orange dashed lines show the slope expected for the analytic scaling in Eq. (12) (slope = 2) and the approximate scaling observed here (slope = 1) for finite-duration pulses. Thick horizontal lines indicate δ=1 and 10 km; gray shading represents δ≤10 km, e.g., the lowermost 10 km of the CLM.


4 Discussion

The model scenarios considered above are idealized and therefore limited in their representation of the complexities of deformation and fluid–rock interactions within the CLM. In particular, the effective thermal properties and the geometry of the fluid- and melt-rich channels are abstracted into a single number, the heat transfer coefficient, K, strongly controlled by the channel spacing, d. Sinuosity and other aspects of the geometry of channelization are abstracted, and the details of processes at and below the scale of an average channel spacing, d, are ignored. Instead, the focus here is on the effective behavior at mesoscopic spatial scales d. Even at these scales, this formulation ignores spatial variations in transport, including variability in the channel volume fraction ϕ, in-channel velocity vchannel, and effective heat transfer coefficient K. Time-dependent variability in transport is also ignored, e.g., feedbacks due to possible phase changes during disequilibrium heating or cooling which would affect the geometry of the channels (Keller and Suckale2019). Finally, this 1D model ignores the 3D nature of relative motion between material inside and outside channel walls even on the mesoscale (d).

Given these limitations, the models above are a way to frame first-order questions and develop arguments related to the consequences of disequilibrium heating, particularly when the behavior is dominated by downstream effects in the direction of transport. Taking the model domain to be analogous to the lowermost lithosphere, where melt or fluid transport may be channelized (Fig. 8), x<0 corresponds to a melt-rich sub-lithospheric region (e.g., a decompaction layer, Holtzman and Kendall2010), whereas the domain x>0 represents an initially sub-solidus lowermost CLM, and x=0 is the initial LAB (Fig. 8). Melt infiltration into the lithosphere may be episodic, controlled by timescales associated with transport from the melt-generation zone to the LAB (e.g., Scott and Stevenson1984; Wiggins and Spiegelman1995), processes of fracturing and crystallization in a dike boundary layer (e.g., Havlin et al.2013) and melt supply from a deeper region of melt production (e.g., Lamb et al.2017).

Figure 8Implications for a thermal re-working zone (TRZ) that forms a modified layer at the lowermost CLM as a result of disequilibrium heating. (a) Illustration of a specific set of parameters that lead to TRZ growth to a steady-state width of δ≈10 km, after 0.1 to 1 Myr: duration of melt infiltration τ=101 to 102 kyr, channel spacing d=500 to 1000 , and in-channel material velocity vchannel=1 m yr−1. The TRZ is characterized by an upward-decreasing degree of disequilibrium (indicated by the color). Before reaching its final width, the TRZ transiently grows at a relatively fast rate, Vdiseqm>10 km Myr−1. (b) Illustration of an interpretation of spatially variable TRZ width, δ(x), as due to lateral juxtaposition of multiple regions of dominantly vertical (1D) transport, but with spatially variable d, ϕ, vchannel, and τp. In this case, δ1>δ2, which could be due to d1<d2, or ϕ1>ϕ2 or v1>v2, or τp1>τp2, or some combination of these.


Three key results emerge from the models above: (1) disequilibrium heating, estimated using the heat transfer coefficient, may be a significant portion of the heat budget at the LAB and the lowermost CLM; (2) there is a material-dependent velocity associated with transient disequilibrium heating; and (3) there is a region of spatiotemporally varying disequilibrium heat exchange, a thermal reworking zone (TRZ). Below I discuss each of these within the context of episodic melt infiltration into the CLM in an intra-plate setting, specifically the Basin and Range province of the western US where deformation and 3D melt transport may be simplified by neglecting plate-boundary effects. In this case, dominantly vertical heat transport within a slowly deforming lithosphere is a reasonable first-order assumption.

i. Disequilibrium heating and the heat budget at the LAB. The relative importance of disequilibrium heat exchange at the LAB may be established by considering the effective heat transfer coefficient, K, and the parameter which most strongly controls it, namely the average spacing of channels, d. For the material parameters in Table 1, and channel volume fraction ϕ≈0.1, channel spacing of d=500 to 1000 m, K is of order 10−5 W m−3 K−1 (Sect. 2.1), corresponding to Péclet numbers on the order of 101. Physically, K corresponds to across-channel heat transfer per unit time, per unit volume, and per unit difference in temperature (Schumann1929). Therefore, if we assume that the average thermal contrast in the TRZ is roughly 2 % to 5 % of ΔT (e.g., (Tf-Ts)max in Fig. 6), for ΔT=100 K excess temperature of the infiltrating melt, disequilibrium heat exchange might contribute around 10−4 to 10−5 W m−3 to the heat budget at the LAB. This is a conservative estimate, given that the temperature difference between magma and the surrounding material may be larger (e.g., in Lherz the inferred contrast is >200 K, (Soustelle et al.2009), and up to 1000 K in crust, Lesher and Spera2015). Similarly, plume excess temperatures are estimated to be as large as 250 K (Wang et al.2015).

To put this in perspective, we now compare this estimated heat budget to the heat budget due to deposition of latent heat during crystallization of melt transported in channels in the lithosphere (e.g., Havlin et al.2013; ReesJones et al.2018). The contribution from freezing of melt may be estimated using scaling arguments made in Havlin et al. (2013). Assuming that melt and rock are in equilibrium, Havlin et al. (2013) estimate that the heat released by a crystallization front would contribute around ρHSdike, where ρ is the melt density, H is the latent heat of crystallization, and Sdike is a volumetric flow rate out of a decompacting melt-rich LAB boundary layer due to diking. For a representative dike porosity of ϕ=0.1 within the dike, Havlin et al. (2013) estimate Sdike≈0.2 mm3 s−1. Taking ρ=3000 kg m3, and H=3×105 J kg−1, the heat source due to the moving crystallization front would be around 102 W per dike. If we assume that this heating takes place within a volume that is roughly the dike height × dike spacing × dike length, we can determine the volumetric power generated due to crystallization. For example, assuming dike heights of about 1 km and dike spacing large enough for non-interacting dikes (as estimated by Havlin et al.2013, a porosity of 0.1 would require a dike spacing of ∼1 km), the heat source due to a crystallizing dike boundary layer would be <10-4 W m−3 (per unit length along strike). These arguments corroborate the idea that disequilibrium heat exchange during melt–rock interaction could be an important portion of the heat budget at the LAB as compared to other expected processes, such as heating due to crystallization of melt in channels (see also ReesJones et al.2018).

ii. Progression of a disequilibrium heating zone or front at a rate Vdiseqm. The disequilibrium heating front is associated with a migration rate Vdiseqm that is less than the in-channel material velocity, vchannel (Eq. 11). Therefore, Vdiseqm limits the rate at which the lowermost CLM may be modified by thermal disequilibrium during migration of a disequilibrium front (e.g., Figs. 5b, 6). Although in this work I considered a fixed vchannel=1 m yr−1, note that Vdiseqm (Eq. 11) is independent of temperature contrast between the CLM and infiltrating melt and, for fixed material properties and channel volume fraction ϕ, depends linearly on vchannel. However, Eq. (11) also illustrates a tradeoff between vchannel and ϕ. Assuming a channel volume fraction of ϕ=0.1 at the LAB, and material properties in Table 1, Vdiseqm would be around 10 % of the in-channel velocity (Fig. 4c). For in-channel velocity in the range of 0.01 to 1 m yr−1, the disequilibrium heating front at the LAB would migrate upward at a rate of ≈1 to 102 km Myr−1 (Fig. 8) This is comparable to rates of CLM thinning predicted by heating due to the upward motion of a dike boundary layer (1 to 6 km Myr−1 in Havlin et al.2013). Interestingly, an upward-moving disequilibrium heating zone with Vdiseqm≈1 to 102 km Myr−1 also brackets the 10–20 km Myr−1 rate of upward migration of the LAB inferred from the pressure and temperature of last equilibration of Cenozoic basalts in the Big Pine volcanic field in the western US (Plank and Forsyth2016). An implication of the models here, therefore, is that disequilibrium heating may produce lithosphere modification at geologically relevant spatial and temporal scales provided that the material velocity in channels at the LAB is on the order of 10−1 to 1 m yr−1 (e.g., Fig. 8Rutherford2008); higher transport rates would require lower ϕ to drive a similar rate of CLM modification.

iii. Thermal reworking zone (TRZ). A key result that may be relevant to the evolution of the LAB is that episodic infiltration of melts that are hotter than the surrounding CLM would lead to a finite region of disequilibrium heating within a thermal reworking zone or TRZ. The TRZ would undergo a phase of transient widening (at a rate given by Vdiseqm), reaching a maximum width δ that is proportional to ϕvchan(τ/d)n (n≈1 for individual perturbations, but <2 for multiple; Fig. 7; also Fig. 5d). Here d is a characteristic scale of channelization and τ is a timescale associated with the episodicity of melt infiltration. This scaling gives us a way to conceptualize the modification of the lowermost CLM as a zone that may encompass a variable thickness TRZ, depending on variability in transport velocity and in the timescale of melt infiltration (Fig. 8b). Regions where the timescale of episodic melt infiltration is longer are predicted to have a thicker zone of modification at the LAB (Fig. 7). For example, for ϕ≈0.12, a channel spacing of d=500 m, disequilibrium heating by melt pulses that last around 10 kyr implies a maximum thickness of roughly 10 km for the zone of modification (Fig. 7; also Fig. 5d). In this scenario, the TRZ grows to this maximum width over a timescale governed by δ/Vdiseqm; for Vdiseqm=10 km Myr−1, which corresponds to melt velocity of roughly 0.1 m yr−1 (see (ii) above), the 10 km wide TRZ would be established within about 1 Myr (Fig. 8), comparable to rates of CLM modification inferred from observations in Plank and Forsyth (2016).

These scaling arguments lead to the idea that perhaps the TRZ represents a zone of thermal modification at the base of the CLM that may also correspond to (or encompass) a zone of rheologic weakening and/or in situ melting if the infiltrating fluids are hotter than the ambient material. The dynamic evolution of the LAB during episodic pulses of melt infiltration is beyond the scope of the simple models above (which assume a stationary, undeforming matrix). However, assuming mantle material obeys a temperature and pressure-dependent viscosity scaling relation such as in Hirth and Kohlstedt (2003) at an LAB depth of about 75 km, where the ambient mantle is cooler than the dry solidus (e.g., 1100 C + 3.5 C km−1;, Plank and Forsyth2016, and assuming a wet dislocation creep mechanism), we would expect a viscosity reduction by a factor 1/1.4 to 1/12 during heating (e.g., 1/2.3 for a temperature increase of 20 K (0.2ΔT for a 100 K perturbation amplitude)). This effect is weaker, but still important for a deeper LAB; e.g, at 125 km depth, the viscosity reduction would be a factor 1/1.2 for a temperature increase of 20 K.

Geochemical evidence from Cenozoic basalts from the western US, particularly space–time variations in volcanic rock Ta/Th and Nd isotopic compositions, suggests that the timescale of modification and removal of the lowermost CLM is on the order of 101 Myr (Farmer et al.2020). These authors argue that the observed transition from low to intermediate to high Ta/Th ratios indicates a change from arc- or subduction-related magmatism, to magmatism associated with in situ melting of a metasomatized CLM (the “ignimbrite flare-up”), to magmatism due to decompression and upwelling after removal of the lowermost CLM. At a minimum, the observed timescale of the transition in Ta/Th ratios in volcanic rocks (101 Myr) in the western US should be comparable to the timescales of degradation of the CLM. If correct, these interpretations and observations are promising and provide an important avenue for exploring the role of thermal and chemical disequilibrium during melt–rock interaction and destabilization of the CLM in an intra-plate setting.

5 Conclusions

In summary, I have presented arguments supporting the role of disequilibrium heating in the modification of the base of the continental lithospheric mantle (CLM) during melt infiltration into and across the lithosphere–asthenosphere boundary (LAB). Infiltration of pulses of hotter-than-ambient material into the LAB should establish a thermal reworking zone (TRZ) associated with disequilibrium heat exchange. The spatial and temporal scales associated with the establishment of the TRZ are comparable to those for CLM modification inferred from geochemical and petrologic observation intra-plate settings, e.g., the western US. Disequilibrium heating may contribute around 10−5 W m−3 to the heat budget at the LAB and, for transport velocity of 0.1 to 1 m yr−1 in channels that are roughly 102 m apart, a 10 km wide TRZ may be established within 1 Myr. Disequilibrium heating during melt infiltration may therefore be an important process for modifying the CLM. Further work is needed to explore its role in the rheologic weakening that must precede mobilization (and possibly removal) of the lowermost CLM.

Code availability

The codes (written in MATLAB) used for the non-dimensional system presented are available at (Roy2022).

Data availability

This work does not use any data directly. Instead, published observations are used to support the interpretations here.

Competing interests

The author has declared that there are no competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work benefited from detailed, thoughtful, and constructive reviews from Harro Schmeling and an anonymous referee. The author also wishes to thank topical editor, Juliane Danneberg. I have greatly enjoyed conversations with Lang Farmer, Alisha Clark, and Rick Carlson regarding applications of the simple 1D model presented here to the western US and generally to melt–rock interactions in the real world. I am also grateful for discussions with Chris Havlin and Ben Holtzman during early stages of thinking about the 1D Schumann model. Funding sources for developing the model and ideas presented include NSF EAR-0952325, EAR-2120812, an Aspen Center for Physics Fellowship, and a Women in STEM award from UNM-Advance. This study draws upon previously published observations, and data and observations were not generated for this research.

Financial support

This research has been supported by the Directorate for Geosciences (grant nos. 2120812 and 0952325) and the Aspen Center for Physics (fellowship to Mousumi Roy grant).

Review statement

This paper was edited by Juliane Dannberg and reviewed by Harro Schmeling and one anonymous referee.


Aharonov, E., Whitehead, J., Kelemen, P., and Spiegelman, M.: Channeling instability of upwelling melt in the mantle, J. Geophys. Res.-Solid Ea., 100, 20433–20450, 1995. a

Bo, T., Katz, R. F., Shorttle, O., and Rudge, J. F.: The melting column as a filter of mantle trace-element heterogeneity, Geochem. Geophy. Geosy., 19, 4694–4721, 2018. a, b

Bodinier, J. L., Vasseur, G., Vernie'res, J., Dupuy, C., and Fabrie's, J.: Mechanisms of mantle metasomatism – geochemical evidence from the Lherz orogenic peridotite, J. Petrol., 31, 597–628, 1990. a

Bodinier, J.-L., Garrido, C. J., Chanefo, I., Brugier, O., and Gervilla, F.: Origin of Pyroxenite-Peridotite Veined Mantle by Refertilization Reactions: Evidence from the Ronda Peridotite (Southern Spain), J. Petrol., 49, 999–1025, 2008. a

Carlson, R. W., Irving, A. J., Schulzec, D. J., and Hearn Jr., B. C.: Timing of Precambrian melt depletion and Phanerozoic refertilization events in the lithospheric mantle of the Wyoming Craton and adjacent Central Plains Orogen, Lithos, 77, 453–472, 2004. a

Chevalier, L. and Schmeling, H.: Thermal non-equilibrium of porous flow in a resting matrix applicable to melt migration: a parametric study, Solid Earth, 13, 1045–1063,, 2022. a, b, c, d, e, f

Dixon, A. G. and Cresswell, D. L.: Theoretical prediction of effective heat transfer parameters in packed beds, AIChE Journal, 25, 663–676, 1979. a, b, c, d, e, f

Dullien, F. A. L.: Porous Media: Fluid Transport and Pore Structure, Academic Press (New York), ISBN 0-12-223651-3, 1979. a, b

England, P. C. and Katz, R. F.: Melting above the anhydrous solidus controls the location of volcanic arcs, Nature, 467, 700–704, 2010. a

Farmer, G. L., Fritz, D., and Glazner, A. F.: Identifying Metasomatized Continental Lithospheric Mantle Involvement in Cenozoic Magmatism from Ta/Th Values, Southwestern North America, Geochem. Geophy. Geosy., 21, e2019GC008499,, 2020. a, b, c

Gao, S., Rudnick, R. L., Carlson, R. W., McDonough, W. F., and Liu, Y.-S.: Re-Os evidence for replacement of ancient mantle lithospherebeneath the North China craton, Earth Planet. Sc. Lett., 198, 307–322, 2002. a

Handley, D. and Heggs, P. J.: Momentum and heat transfer mechanisms in regular shaped packings, Trans. Inst. Chem. Engrs., 46, 251–259, 1968. a, b

Hauri, E. H.: Melt migration and mantle chromatography, 1: simplified theory and conditions for chemical and isotopic decoupling, Earth Planet. Sc. Lett., 153, 1–19, 1997. a, b

Havlin, C., Parmentier, E., and Hirth, G.: Dike propagation driven by melt accumulation at the lithosphere-asthenosphere boundary, Earth Planet. Sc. Lett., 376, 20–28,, 2013. a, b, c, d, e, f, g, h

Hirth, G. and Kohlstedt, D. L.: Rheology of the Upper Mantle and the Mantle Wedge: A View from the Experimentalists, Geophys. Monog. Series, 138, 83–105, 2003. a

Holtzman, B. K. and Kendall, J.-M.: Organized melt, seismic anisotropy, and plate boundary lubrication, Geochem. Geophys. Geosyst., 11, Q0AB06,, 2010. a, b

Hopper, E., Gaherty, J. B., Shillington, D. J., Accardo, N. J., Nyblade, A. A., Holtzman, B. K., Havlin, C., Scholz, C. A., Chindandali, P. R. N., Ferdinand, R. W., Mulibo, G. D., and Mbogoni, G.: Preferential localized thinning of lithospheric mantle in the melt-poor Malawi Rift, Nat. Geosci., 13, 584–589, 2020. a

Keller, T. and Suckale, J.: A continuum model of multi-phase reactive transport in igneous systems, Geophys. J. Int., 219, 185–222, 2019. a, b, c

Kenyon, P. M.: Trace element and isotopic effects arising from magma migration beneath mid-ocean ridges, Earth Planet. Sc. Lett., 101, 367–378, 1990. a, b, c, d

Kuznetsov, A. V.: An investigation of a wave of temperature difference between solid and fluid phases in a porous packed bed, Int. J. Heat Mass Tran., 37, 3030–3033, 1994. a, b, c, d

Kuznetsov, A. V.: Comparison of the waves of temperature difference between the solid and fluid phases in a porous slb and in a semi-infinite porous body, Int. Commun. Heat Mass, 22, 499–506, 1995a. a, b

Kuznetsov, A. V.: An Analytical Solution for Heating a Two-Dimensional Porous-Packed Bed by a Non-Thermal Equilibrium Fluid Flow, Appl. Sci. Res., 55, 83–93, 1995b. a, b

Kuznetsov, A. V.: Analysis of heating a three-dimensional porous bed utilizing the two energy equation model, Heat Mass Transfer, 31, 173–178, 1996. a, b

Lamb, S., Moore, J. D. P., Smith, E., and Stern, T.: Episodic kinematics in continental rifts modulated by changes in mantle melt fraction, Nature, 547, 84–88, 2017. a

Lenoir, X., Garrido, C., Bodinier, J., Dautria, J., and Gervilla, F.: The recrystallization front of the ronda peridotite: Evidence for melting and thermal erosion of subcontinental lithospheric mantle beneath the alboran basin, J. Petrol., 42, 141–158, 2001. a

LeRoux, V., Bodinier, J., Tommasi, A., Alard, O., Dautria, J., Vauchez, A., and Riches, A.: The lherz spinel lherzolite: Refertilized rather than pristine mantle, Earth Planet. Sc. Lett., 259, 599–612, 2007. a, b, c, d, e

LeRoux, V., Bodinier, J.-L., Alard, O., O'Reilly, S. Y., and Griffin, W. L.: Isotopic decoupling during porous melt flow: A case-study in the Lherz peridotite, Earth Planet. Sc. Lett., 279, 76–85, 2008. a, b, c, d, e

Lesher, C. E. and Spera, F. J.: Encylcopedia of Volcanoes, 2nd Edn., Chapter 5 – Thermodynamic and Transport Properties of Silicate Melts and Magma, 113–141, Elsevier,, 2015. a, b, c, d, e, f, g

Menzies, M., Xu, Y., Zhang, H., and Fan, W.: Integration of geologygeophysics and geochemistry: A key to understanding the North China Craton, Lithos, 96, 1–21, 2007. a

O'Reilly, S. Y., Griffin, W. L., Djomani, Y. H. M., and Morgan, P.: Are lithospheres forever? Tracking changes in the subcontinental lithospheric mantle through time, GSA Today, 11, 4–10, 2001. a

Pec, M., Holtzman, B. K., Zimmerman, M. E., and Kohlstedt, D. L.: Reaction Infiltration Instabilities in Mantle Rocks: an Experimental Investigation, J. Petrol., 58, 979–1004, 2017. a

Perrin, A., Goes, S., Prytulak, J., Davies, D. R., Wilson, C., and Kramer, S.: Reconciling mantle wedgethermal structure with arc lavathermobarometric determinations inoceanic subduction zones, Geochem. Geophy. Geosy., 17, 4105–4127, 2016. a, b

Plank, T. and Forsyth, D. W.: Thermal structure and melting conditions in the mantle beneath the Basin and Range province from seismology and petrology, Geochem. Geophy. Geosy., 17, 1312–1338, 2016. a, b, c, d, e, f, g

ReesJones, D. W., Katz, R. F., Tian, M., and Rudge, J. F.: Thermal impact of magmatism in subduction zones, Earth Planet. Sc. Lett., 481, 73–79, 2018. a, b, c, d, e

Roy, M.: Thermal disequilibrium during melt-transport: Implications for the evolution of the lithosphere-asthenosphere boundary, arXiv [preprint],, 2020. a, b, c, d, e

Roy, M.: mousumiroy-unm/diseqm22: Roy_thermal_diseqm_repo (v1.0.0), Zenodo [code],, 2022. a

Roy, M., Gold, S., Johnson, A., Orozco, R. O., Holtzman, B. K., and Gaherty, J.: Macroscopic coupling of deformation and melt migration at continental interiors, with applications to the Colorado Plateau, J. Geophys. Res., 21, 3762–3781,, 2016. a

Rutherford, M.: Magma ascent rates, Rev. Mineral. Geochem., 69, 241–271, 2008. a, b

Saal, A. E., Takazawa, E., Frey, F. A., Shimizu, N., and Hart, S. R.: Re-Os isotopes in the Horoman peridotite: Evidence for refertilization?, J. Petrol., 42, 25–37, 2001. a

Sakamaki, T., Suzuki, A., Ohtani, E., Terasaki, H., Urakawa, S., Katayama, Y., ichi Funakoshi, K., Wang, Y., Hernlund, J. W., and Ballmer, M. D.: Ponded melt at the boundary between the lithosphere and asthenosphere, Nat. Geosci., 6, 1041–1044,, 2013. a

Schmeling, H., Marquart, G., and Grebe, M.: A porous flow approach to model thermal non-equilibrium applicable to melt migration, Geophys. J. Int., 212, 119–138, 2018. a, b, c, d, e, f, g, h, i

Schumann, T. E.: Heat transfer: a liquid flowing through a porous prism, J. Franklin Inst., 208, 405–416, 1929. a, b, c, d, e, f, g, h, i, j, k

Scott, D. R. and Stevenson, D. J.: Magma solitons, Geophys. Res. Lett., 11, 1161–1164, 1984.  a

Soustelle, V., Tommasi, A., Bodinier, J. L., Garrido, C. J., and Vauchez, A.: Deformation and Reactive MeltTransport in the Mantle Lithosphere above a Large-scale Partial Melting Domain: the Ronda Peridotite Massif, Southern Spain, J. Petrol., 50, 1235–1266, 2009. a, b, c

Spiga, G. and Spiga, M.: A rigorous solution to a heat transfer two-phase model in porous media and packed beds, Int. J. Heat Mass, 24, 355–364, 1981. a, b, c, d, e

Vauchez, A. and Garrido, C. J.: Seismic properties of an asthenospherized lithospheric mantle: constraints from lattice preferred orientations in peridotite from the Ronda massif, Earth Planet. Sc. Lett., 192, 235–249, 2001. a

Wallner, H. and Schmeling, H.: Numerical models of mantle lithosphere weakening, erosion and delamination induced by melt extraction and emplacement, Int. J. Earth Sci. (Geologische Rundschau), 105, 1741–1760, 2016. a, b

Wang, H., van Hunen, J., and Pearson, D. G.: The thinning of subcontinental lithosphere: The roles of plume impact and metasomatic weakening, Geochem. Geophy. Geosy., 16, 1156–1171,, 2015. a, b

Wenker, S. and Beaumont, C.: Can metasomatic weakening result in the rifting of cratons?, Tectonophysics, 746, 3–21, 2017. a

Wiggins, C. and Spiegelman, M.: Magma migration and magmatic solitary waves in 3-D, Geophys. Res. Lett., 22, 1289–1292, 1995. a

Short summary
This study investigates one of the key processes that may lead to the destruction and destabilization of continental tectonic plates: the infiltration of buoyant, hot, molten rock (magma) into the base of the plate. Using simple calculations, I suggest that heating during melt–rock interaction may thermally perturb the tectonic plate, weakening it and potentially allowing it to be reshaped from beneath. Geochemical, petrologic, and geologic observations are used to guide model parameters.