**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

Mousumi Roy

**Mousumi Roy**Mousumi Roy

- Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87106, USA

- Department of Physics and Astronomy, University of New Mexico, Albuquerque, NM 87106, USA

**Correspondence**: Mousumi Roy (mroy@unm.edu)

**Correspondence**: Mousumi Roy (mroy@unm.edu)

Received: 31 Jan 2022 – Discussion started: 17 Feb 2022 – Revised: 28 Jul 2022 – Accepted: 09 Aug 2022 – Published: 05 Sep 2022

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 (*v*_{channel}) 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 $\left[\mathit{\varphi}{v}_{\mathrm{channel}}(\mathit{\tau}/d{)}^{n}\right]$, 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 *v*_{channel}=1 m yr^{−1}, channel spacings *d*≈10^{2} m, and timescales of episodic melt infiltration *τ*≈10^{1} kyr, the steady-state width of the TRZ in the lowermost CLM is ≈10 km.
(3) Within the TRZ, disequilibrium heat exchange may contribute $\approx {\mathrm{10}}^{-\mathrm{5}}$ W m^{−3} to the LAB heat budget.

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 Beaumont, 2017; Plank and Forsyth, 2016; 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 Katz, 2010; 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 Garrido, 2001; 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 ($\approx \mathrm{230}\phantom{\rule{0.125em}{0ex}}{}^{\circ}$C 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 Forsyth, 2016), 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 Forsyth, 2016), 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 Kendall, 2010) 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 Schmeling, 2022). 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 Schmeling, 2022; Keller and Suckale, 2019; Wallner and Schmeling, 2016; 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 Schmeling, 2022) in order to explore thermal equilibration over long timescales of melt infiltration ($>{\mathrm{10}}^{\mathrm{3}-\mathrm{4}}$ years) into the lowermost 1–10 km of the CLM (e.g., stage 3, large Péclet numbers in Chevalier and Schmeling, 2022).

The 1D model abstracts the complex geometry of the melt–rock interface and therefore differs from other descriptions of disequilibrium heat exchange (Wallner and Schmeling, 2016; Schmeling et al., 2018; Keller and Suckale, 2019). Similar to chemical transport models that use a linear driving term for chemical disequilibrium (e.g., Hauri, 1997; Kenyon, 1990; Bo et al., 2018), the models below assume a linear thermal driving term (e.g., Schumann, 1929; Kuznetsov, 1994; Spiga and Spiga, 1981). In this study, however, only the role of thermal disequilibrium is included and I ignore important processes such as chemical disequilibrium (e.g., Hauri, 1997; Kenyon, 1990; 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 Forsyth, 2016; Farmer et al., 2020).

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, *v*_{channel}, 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 Roy, 2020).
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, *T*_{s}, and within channels, *T*_{f}:

where *v*_{channel} 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, *c*_{f} and *c*_{s} are the heat capacities per unit volume (heat capacitances) at constant pressure (*c*_{f}=*c*_{pf}*ρ*_{f} and *c*_{s}=*c*_{ps}*ρ*_{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 Spiga, 1981; Kuznetsov, 1994, 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 *v*_{channel} 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.

In this model, the terms involving the thermal contrast (*T*_{f}−*T*_{s}) 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, ${t}_{f}=\mathrm{1}/{K}_{f}=\mathit{\varphi}{c}_{f}/K$, and outside channels, ${t}_{s}=\mathrm{1}/{K}_{s}=(\mathrm{1}-\mathit{\varphi}){c}_{s}/K$.
Additionally, a characteristic length scale emerges out of the relative motion across channel walls, ${v}_{\mathrm{channel}}/{K}_{f}$.
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, ${K}_{s}=K/\left({c}_{s}\right(\mathrm{1}-\mathit{\varphi}\left)\right)$ and ${K}_{f}=K/\left({c}_{f}\mathit{\varphi}\right)$ in Eqs. (1) and (2).
Note that *K*_{f} and *K*_{s} 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 Schumann, 1929, 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), *a*_{sf}, 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 ${S}_{\mathrm{0}}=\mathrm{6}/p$, so ${a}_{\mathrm{sf}}={S}_{\mathrm{0}}(\mathrm{1}-\mathit{\varphi})=\mathrm{6}(\mathrm{1}-\mathit{\varphi})/p$ (Dullien, 1979).
This sets a limit for channels with channel spacing *d*, where we shall assume that the specific surface area is ${a}_{\mathrm{sf}}\approx A(\mathrm{1}-\mathit{\varphi})/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 Cresswell, 1979; Schmeling et al., 2018).
Whereas the specific wall area is a geometric factor, the effective conductivity of the medium depends on the Nusselt number, *N**u*.
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):

where *β*=10 for spherical matrix grains, 8 for cylinders, and 6 for slabs (Dixon and Cresswell, 1979).
Therefore, the range of *β*=6 to 10 represents the highly channelized vs. porous flow end-members.
For slow flows (Reynolds number *R**e*≪100), Handley and Heggs (1968) argue that *N**u* ranges from 0.1 to 12.4 (Dixon and Cresswell, 1979) (Table 1).
The relevant quantity that determines *K* is an effective “conductance” ${C}_{\mathrm{eff}}/d$, so that $K={C}_{\mathrm{eff}}{a}_{\mathrm{sf}}/d$,

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 Heggs, 1968)Dixon and Cresswell (1979)Dullien (1979)LeRoux et al. (2008)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, ${\mathit{\lambda}}_{f}={\mathit{\phi}}_{\mathrm{f}}{\mathit{\lambda}}_{\mathrm{melt}}+(\mathrm{1}-{\mathit{\phi}}_{\mathrm{f}}){\mathit{\lambda}}_{\mathrm{grain}}$, where *λ*_{melt} and *λ*_{grain} are the values for, say, basaltic melt and peridotitic grains.
Similarly, outside channels the volume average would be ${\mathit{\lambda}}_{s}={\mathit{\phi}}_{\mathrm{s}}{\mathit{\lambda}}_{\mathrm{melt}}+(\mathrm{1}-{\mathit{\phi}}_{\mathrm{s}}){\mathit{\lambda}}_{\mathrm{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 Spera, 2015), *λ*_{grain}=2.5 W m^{−1} K^{−1}) for the solid grains, and taking *N**u*=0.1 to 12.4, *β*=6, 8, or 10, and *A*=6, we find that the effective conductivity *C*_{eff}, *a*_{sf}, 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), $K\sim {d}^{-\mathrm{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 Schmeling, 2022).
In the models in Sect. 3 I set *β*=6 and *A*=2 as suggested for 1D channels (Dixon and Cresswell, 1979; 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, *K*∼*C*_{eff} $/$ (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 $K\sim {d}^{-\mathrm{2}}$, our models must be confined to timescales that are large relative to the material-dependent thermal response timescale (e.g., $\mathrm{1}/{K}_{s}$), a requirement that is met in all thermal perturbations considered 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*=10^{0} to 10^{3} m (1 m to 1 km channel spacing).
The corresponding range of the heat transfer coefficient in the models is therefore $K\approx {\mathrm{10}}^{-\mathrm{5}}$ (large *d*) to 10^{1} (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*>10^{2} m in order to ensure an effective Péclet number *P**e*>10.
Taking typical parameters *ϕ*=0.1, *d*=1000 m, and corresponding *K* values in Fig. 2, typical timescales of response are ${t}_{s}=\mathrm{1}/{K}_{s}\approx \mathrm{9800}$ years and ${t}_{f}=\mathrm{1}/{K}_{f}\approx \mathrm{1200}$ 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 Spiga, 1981), we define the normalized relative temperature, ${T}_{f}^{\prime}=({T}_{f}-{T}_{\mathrm{0}})/\mathrm{\Delta}T$ and ${T}_{s}^{\prime}=({T}_{s}-{T}_{\mathrm{0}})/\mathrm{\Delta}T$, where *T*_{0} is reference temperature and Δ*T* is a temperature perturbation (described below).
We also introduce the dimensionless position, ${x}^{\prime}=x{K}_{f}/{v}_{\mathrm{channel}}$, a dimensionless time, ${t}^{\prime}={K}_{s}t$, and the weighted heat capacitance ratio, $\mathit{\zeta}={K}_{s}/{K}_{f}=\mathit{\varphi}{c}_{f}/(\mathrm{1}-\mathit{\varphi}){c}_{s}$ (from Eq. (5), we see that $\mathrm{1}/\mathit{\zeta}$ 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 ${D}_{f}={D}_{s}=\mathrm{0}$ (Spiga and Spiga, 1981; Kuznetsov, 1994, 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, *D*_{f} and *D*_{s}.
Using the definition of *ζ*, we can further simplify the coefficient of the diffusion term *D*_{s} and express it in terms of *D*_{f} as

It is clear that for a given temperature difference $({T}_{f}^{\prime}-{T}_{s}^{\prime})$ the behavior of Eqs. (5) and (6) is governed by *D*_{f} and *ζ* ($\mathrm{1}/\mathit{\zeta}$ 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 *v*_{channel}, the characteristic length scale ${v}_{\mathrm{channel}}/{K}_{f}={v}_{\mathrm{channel}}\left(\mathit{\varphi}{c}_{f}\right)/K$, divided by the thermal diffusivity of the channel material ${\mathit{\lambda}}_{f}/{c}_{f}$,

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

Alternatively, one may also define an effective Péclet number as the product of ${v}_{\mathrm{channel}}^{\mathrm{2}}$, a characteristic timescale $\mathrm{1}/{K}_{s}={c}_{s}(\mathrm{1}-\mathit{\varphi})/K$, divided by the thermal diffusivity of material outside channels ${\mathit{\lambda}}_{s}/{c}_{s}$,

These definitions are limited to consideration of channelized flow and more complex dependence between *P**e*, and *ϕ* is suggested in the case of tubes or pores (e.g., Chevalier and Schmeling, 2022).
Using values of *N**u*=12, *β*=6, *A*=2, with *ϕ* and material properties in Table 1, for a given *v*_{channel}, *d*, *P**e*_{1} is ≈2 to 3 times larger than *P**e*_{2}.
Hereafter, I use the definition of *P**e*_{2} when referring to Péclet number *P**e*.

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, $P{e}_{\mathrm{channel}}={v}_{\mathrm{channel}}d{c}_{f}/{\mathit{\lambda}}_{f}$.
By using the characteristic times $\mathrm{1}/{K}_{f}$ or $\mathrm{1}/{K}_{s}$, 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}^{\prime}=[\mathrm{0},L]$ with *L*=10^{2} to 10^{3} is discretized typically with *N*=1000–5000 elements, and the maximum time for each run is ${t}_{\mathrm{max}}^{\prime}>L\mathit{\zeta}$ (for thermal perturbations considered below that have a characteristic time *τ*, I use ${t}_{\mathrm{max}}^{\prime}>L\mathit{\zeta}+\mathrm{2}{K}_{s}\mathit{\tau}$).
For a given element size d${x}^{\prime}=L/N$, second-order finite differences are used for all spatial derivatives.
The time step d*t*^{′} is chosen (empirically) to be small enough to avoid numerical dispersion (for *P**e*>1, it is sufficient to choose $\mathrm{d}{t}^{\prime}=\mathrm{d}{x}^{\prime}\times {\mathrm{10}}^{-\mathrm{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 (${D}_{f}={D}_{S}=\mathrm{0}$) were benchmarked (using leap-frog with upwind differencing for the spatial derivatives) against analytic solutions in Schumann (1929).

In the models below, I consider a 1D domain (*x*≥0) that is initially in equilibrium ${T}_{s}^{\prime}={T}_{f}^{\prime}=\mathrm{0}$; material outside the channels is stationary, and only the in-channel material is moving, at speed *v*_{channel}.
At *x*=0 and $t={\mathrm{0}}^{+}$, the in-channel material is subjected to a thermal perturbation in *T*_{f}.
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 *P**e*>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 *T*_{f} disturbs the initial steady-state condition in the domain starting at $t={\mathrm{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 *x*_{pert}=*v*_{channel}*t* (or ${x}_{\text{pert}}^{\prime}={t}^{\prime}/\mathit{\zeta}$; where the dimensionless velocity inside channels is $\mathrm{1}/\mathit{\zeta}$; see Sect. 2.2).
(Note that the location of maximum disequilibrium heat exchange $({T}_{f}^{\prime}-{T}_{s}^{\prime}{)}_{\text{max}}$) moves at a different speed, *V*_{diseqm}, 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 *c*_{f}, and *c*_{s}, and the speed of the in-channel material *v*_{channel} (see values in Table 1).
Since the material properties *c*_{s} and *c*_{f} 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 ${T}_{s}={T}_{f}={T}_{\mathrm{0}}$ (${T}_{s}^{\prime}={T}_{f}^{\prime}=\mathrm{0}$), and at *t*=0 the temperature of the in-channel material entering at the inflow *x*=0 is perturbed so that ${T}_{f}(x=\mathrm{0},t>\mathrm{0})={T}_{\mathrm{0}}+\mathrm{\Delta}T$ (or ${T}_{f}^{\prime}=\mathrm{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, ${T}_{s}^{\prime}={T}_{f}^{\prime}=\mathrm{0}$.
Behind this zone, the channels are in equilibrium with the surroundings at the inlet temperature, ${T}_{s}^{\prime}={T}_{f}^{\prime}=\mathrm{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 *V*_{diseqm}, 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 Roy, 2020).
At large *P**e* this disequilibrium heat exchange dominates over axial conduction.

Following an initial lag time (when the maximum thermal contrast is at *x*=0), the disequilibrium zone (marked by the peak in the ${T}_{f}^{\prime}-{T}_{s}^{\prime}$ function, $({T}_{f}^{\prime}-{T}_{s}^{\prime}{)}_{\text{max}}$ in Fig. 4b) migrates inward migration at a steady speed *V*_{diseqm}, a fixed fraction of *v*_{channel} that depends on *ϕ*, the channel volume fraction (Fig. 4a and b).
The ratio ${V}_{\mathrm{diseqm}}/{v}_{\mathrm{channel}}$ varies linearly with *ϕ* (Fig. 4c).
In the near-equilibrium limit (${T}_{f}^{\prime}-{T}_{s}^{\prime}\approx \mathrm{0}$), Kuznetsov (1994) shows that the shape of the temperature difference function (Fig. 4b) approaches a Gaussian function with width that depends on $\sqrt{{t}^{\prime}}$ and the zone of disequilibrium migrates at speed ${v}_{\mathrm{channel}}{c}_{f}/(\mathit{\varphi}{c}_{f}+(\mathrm{1}-\mathit{\varphi}\left){c}_{s}\right)$.
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 $({T}_{f}^{\prime}-{T}_{s}^{\prime}{)}_{\text{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, *v*_{channel},

Although *V*_{diseqm} 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 *v*_{channel}=1 m yr^{−1}, and *d*=500 to 1000 m, which correspond to $K=\mathrm{4.7}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{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.

## 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 $\mathrm{1}/{K}_{s}$, the longest response timescale in the domain (on the order of 10^{3} 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 10^{4} to 10^{6} 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 $\mathrm{1}/{K}_{s}$ penetrate farther into the domain than shorter period oscillations (Fig. 5).
The non-dimensional period, $\mathit{\tau}{K}_{s}=\mathit{\tau}K/(\mathrm{1}-\mathit{\varphi}){c}_{s}$, controls the length scale, *δ*, over which thermal oscillations penetrate into the domain.
The wavelength of these temperature oscillations is set by the period *τ*, *λ*=*v*_{channel}*τ*.
The penetration distance of the oscillations, *δ*, is the maximum width of the TRZ.
At short timescales, when $t\ll \mathit{\delta}/{V}_{\mathrm{diseqm}}$, 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 *V*_{diseqm}, eventually reaching a maximum width, *δ* (akin to a thermal “skin depth”), at time $\mathit{\delta}/{V}_{\mathrm{diseqm}}$.
When $t\gg \mathit{\delta}/{V}_{\mathrm{diseqm}}t$, 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 *V*_{diseqm} (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).

As we might expect, the maximum width of the TRZ, *δ*, is set controlled both by the non-dimensional oscillation period, $\mathit{\tau}/{t}_{s}$ and by the heat transfer coefficient, *K*.
In the absence of the axial diffusion terms, previous work (Spiga and Spiga, 1981; Kenyon, 1990) has shown that

noting that $K\sim {d}^{-\mathrm{2}}$ (2.1) and ${t}_{s}=\mathrm{1}/{K}_{s}={c}_{s}(\mathrm{1}-\mathit{\varphi})/K$; the expression above suggests that, for fixed *v*_{channel}, *δ* is proportional to $(\mathit{\tau}/d{)}^{\mathrm{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 $(\mathrm{1}/\mathrm{2})\left(\mathrm{tanh}\left(\frac{(t-{\mathit{\tau}}_{\mathrm{p}}/\mathrm{2})}{w}\right)-\mathrm{tanh}\left(\frac{(t-\mathrm{3}{\mathit{\tau}}_{\mathrm{p}}/\mathrm{2})}{w}\right)\right)$ 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 ${T}_{f}^{\prime}$ 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 $({T}_{f}^{\prime}-{T}_{s}^{\prime}{)}_{\text{max}}$) migrates into the domain at a rate given by Eq. (11) and decays during transport (Fig. 6b).

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 *v*_{channel}=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, $K\sim {d}^{-\mathrm{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 10^{1} to 10^{2} kyr (Fig. 7).

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 *v*_{channel}, 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 Suckale, 2019).
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 Kendall, 2010), 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 Stevenson, 1984; Wiggins and Spiegelman, 1995), 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).

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 10^{1}.
Physically, *K* corresponds to across-channel heat transfer per unit time, per unit volume, and per unit difference in temperature (Schumann, 1929).
Therefore, if we assume that the average thermal contrast in the TRZ is roughly 2 % to 5 % of Δ*T* (e.g., $({T}_{f}^{\prime}-{T}_{s}^{\prime}{)}_{\text{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 Spera, 2015).
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 *ρ**H**S*_{dike}, where *ρ* is the melt density, *H* is the latent heat of crystallization, and *S*_{dike} 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 *S*_{dike}≈0.2 mm^{3} s^{−1}.
Taking *ρ*=3000 kg m^{3}, and $H=\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ J kg^{−1}, the heat source due to the moving crystallization front would be around 10^{2} 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 $<{\mathrm{10}}^{-\mathrm{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* *V*_{diseqm}.
The disequilibrium heating front is associated with a migration rate *V*_{diseqm} that is less than the in-channel material velocity, *v*_{channel} (Eq. 11).
Therefore, *V*_{diseqm} 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 *v*_{channel}=1 m yr^{−1}, note that *V*_{diseqm} (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 *v*_{channel}.
However, Eq. (11) also illustrates a tradeoff between *v*_{channel} and *ϕ*.
Assuming a channel volume fraction of *ϕ*=0.1 at the LAB, and material properties in Table 1, *V*_{diseqm} 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 10^{2} 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 *V*_{diseqm}≈1 to 10^{2} 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 Forsyth, 2016).
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. 8; Rutherford, 2008); 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 *V*_{diseqm}), reaching a maximum width *δ* that is proportional to $\left[\mathit{\varphi}{v}_{\mathrm{chan}}(\mathit{\tau}/d{)}^{n}\right]$ (*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 $\mathit{\delta}/{V}_{\mathrm{diseqm}}$; for *V*_{diseqm}=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 Forsyth, 2016, and assuming a wet dislocation creep mechanism), we would expect a viscosity reduction by a factor $\approx \mathrm{1}/\mathrm{1.4}$ to $\mathrm{1}/\mathrm{12}$ during heating (e.g., $\approx \mathrm{1}/\mathrm{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 $\approx \mathrm{1}/\mathrm{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 $\mathrm{Ta}/\mathrm{Th}$ and Nd isotopic compositions, suggests that the timescale of modification and removal of the lowermost CLM is on the order of 10^{1} Myr (Farmer et al., 2020).
These authors argue that the observed transition from low to intermediate to high $\mathrm{Ta}/\mathrm{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 $\mathrm{Ta}/\mathrm{Th}$ ratios in volcanic rocks (10^{1} 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.

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 10^{2} 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.

The codes (written in MATLAB) used for the non-dimensional system presented are available at https://doi.org/10.5281/zenodo.6981925 (Roy, 2022).

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

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.

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).

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, https://doi.org/10.5194/se-13-1045-2022, 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 $\mathrm{Ta}/\mathrm{Th}$ Values, Southwestern North America, Geochem. Geophy. Geosy., 21, e2019GC008499, https://doi.org/10.1029/2019GC008499, 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, https://doi.org/10.1016/j.epsl.2013.06.010, 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, https://doi.org/10.1029/2010GC003296, 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, https://doi.org/10.1016/B978-0-12-385938-9.00005-5, 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], https://doi.org/10.48550/arXiv.2009.01496, 2020. a, b, c, d, e

Roy, M.: mousumiroy-unm/diseqm22: Roy_thermal_diseqm_repo (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.6981925, 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, https://doi.org/10.1002/2015JB012149, 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, https://doi.org/10.1038/NGEO1982, 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, https://doi.org/10.1002/2015GC005784, 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