Thermal disequilibrium during melt-transport: Implications for the evolution of the lithosphere-asthenosphere boundary
mmanuscript submitted to
Thermal disequilibrium during melt-transport:Implications for the evolution of thelithosphere-asthenosphere boundary
Mousumi Roy Department of Physics and Astronomy, University of New Mexico
Corresponding author: Mousumi Roy, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] S e p anuscript submitted to Abstract
This study explores melt-infiltration and melt-rock interaction as a means of shaping thecontinental lithosphere from beneath. Specifically, we focus upon the role of thermal dis-equilibrium between melt and the surrounding solid as a means for heating and mod-ifying the continental lithospheric mantle. We use simple pore-network models to esti-mate the importance of heating of the ambient mantle by thermal disequilibrium withmelt transported in channels at the base of the lithosphere. The key physical parame-ters that control the behavior are the fluid volume fraction ( φ ), the relative fluid-solidvelocity ( v fluid ), channel spacing ( d ), and the timescale of episodic melt-infiltration ( τ ).Model results and scaling arguments suggest that for geologically-reasonable values ofthese parameters, disequilibrium heating may contribute more than 10 − W/m to theheat budget at the lithosphere-asthenosphere boundary. We find that there exists a ve-locity scale associated with the transient upward motion of a disequilibrium heating front.During episodic melt-infiltration, the models predict the existence of a thermal rework-ing zone (TRZ) associated with spatio-temporally varying disequilibrium heat exchange.The steady-state width of the TRZ is shown to scale as ∼ (cid:2) φv fluid d τ (cid:3) . The spatio-temporal scales associated with the establishment of the TRZ are comparable with thoseassociated with migration of the lithosphere-asthenosphere boundary as inferred fromgeologic observations within continental intra-plate settings, such as the western US. Continental lithosphere may be profoundly altered by melt-infiltration and meta-somatism, manifest in a range of thermal and chemical changes during melt-rock inter-action. Indeed, melt-rock processes have been implicated in controlling the tectonic sta-bility of continents in a number of settings (Wang et al., 2015; Wenker & Beaumont, 2017)and localities such as the North China craton (O’Reilly et al., 2001; Gao et al., 2002; Men-zies et al., 2007), the Wyoming craton (Carlson et al., 2004), the western US (Plank &Forsyth, 2016; Roy et al., 2016), and the Lherz massif (LeRoux et al., 2007, 2008), amongothers. Melt-enhanced destruction of the lowermost continental lithosphere is also in-ferred as a key process in continental rifting, even within currently magma-poor rifts (e.g.,Hopper et al., 2020). Although evidence for the key role of melt-rock interaction in po-tentially destabilizing continental lithosphere has been steadily growing, there are fewconstraints on the processes involved.This study explores melt-infiltration and melt-rock interaction as a means of shap-ing the continental lithosphere from beneath. Specifically, we focus upon the role of ther-mal disequilibrium between melt and the surrounding solid as a means for heating andmodifying the continental lithospheric mantle (CLM). As evidenced in volcanic eruptions,magma may be transported from depth to the surface in thermal and chemical disequi-librium with the surrounding material. The degree of thermal disequilibrium should in-crease with increasing fluid velocity relative to rock, degree of channelization and there-fore, with decreasing depth within the lithosphere (e.g., Cashman et al., 2017; Schmel-ing et al., 2018). Although disequilibrium thermal and chemical interaction during trans-port is a complex phenomenon (Wallner & Schmeling, 2016; Oliveira et al., 2018; Schmel-ing et al., 2018; Bo et al., 2018; Keller & Suckale, 2019), here we simplify the problemto focus on assessing the temporal and spatial scales over which thermal disequilibriumcan play a role in warming and therefore weakening the lowermost portion of the litho-sphere, near the lithosphere-asthenosphere boundary (LAB). Melts generated deeper inthe mantle, ascending to the LAB and interacting with the CLM are expected to be hot-ter than the surrounding material, transported in channels rather than in percolative flow,with transport Peclet numbers > –2–anuscript submitted to and petrologic evidence for the upward migration of the LAB during melt-rock interac-tion in the western US (e.g., Plank & Forsyth, 2016). The model is based on a simple, 1-D theory of heat exchange in a semi-infinite two-phase medium where we assume that fluid transport only occurs along one dimension(Figure 1; e.g., Schumann (1929)). The domain is assumed to be occupied by a constant
Figure 1.
Cartoon of 1D model with parameters such as: specific heat capacities ( c p ) anddensities ( ρ ), fluid velocity v , fluid volume fraction, φ , and the fluid-solid heat transfer coefficient k . The heat transfer coefficient k is a function of the geometry of the fluid-solid interface andscales as d − (Appendix A); large k corresponds to large fluid-solid contact area per unit volume(e.g., small channel spacing, d ) and vice versa. volume fraction, φ , of fluid moving with a constant (average) velocity v relative to thesurrounding stationary solid (with volume fraction 1 − φ ). The model domain may bethought of as co-moving with the reference frame of a constant-velocity solid matrix. Theresults of this model are therefore applicable to physical situations where heat transportin the fluid-flow direction dominates and where the solid motion is steady.We assume that an unchanging fluid-solid interface separates the phases within thedomain. The detailed geometry of this interface is not specified in this 1-D calculationbut, following Schumann (1929) and subsequent pore-network models (Spiga & Spiga,1981; Kuznetsov, 1994, 1995b, 1995a, 1996), is parametrized by a heat transfer coeffi-cient, k (Figure 1). This coefficient is a proxy for the geometry of the inter-phase inter-face, namely the contact area per unit volume, controlled by the spatial scale of chan-nelization, d . As illustrated in (Figure 1), a large value of k may represent efficient fluid-solid heat exchange as in the case of many channels separated by a small distance. Con-versely, a low value of k would represent inefficient exchange, as in the case of a largercharacteristic length scale between the channels.Following Schumann (1929), we will assume that when the fluid and solid are inthermal disequilibrium, the thermal evolution of the system is governed mainly by heatexchange between the two phases across their interfacial surface. This heat exchange isassumed to dominate over thermal dispersion and axial conductive heat fluxes withinthe solid (e.g., within the walls of channels) and within the fluid (within channels). Ad-ditionally, fluid-solid heat exchange is assumed to be linearly proportional to their lo-cal temperature difference. These arguments lead to coupled equations for the temper-ature of the solid, T s , and of the fluid, T f (Schumann, 1929): ∂T f ∂t + v ∂T f ∂x = − kφc f ( T f − T s ) = − k f ( T f − T s ) (1) –3–anuscript submitted to ∂T s ∂t = k (1 − φ ) c s ( T f − T s ) = k s ( T f − T s ) (2)where c f and c s are the heat capacity per unit volume at constant pressure for solid andfluid, so c f = c p fluid ρ f and c s = c p solid ρ s . Two independent factors on the right handsides specify the timescales of heat exchange in the fluid, 1 /k f = φc f /k , and in the solid,1 /k s = (1 − φ ) c s /k . Similarly, a characteristic length scale emerges out of the fluid mo-tion, v/k f . These characteristic length- and timescales are combinations of five user-specifiedquantities that completely specify the model: the heat transfer coefficient, k , fluid vol-ume fraction φ , the heat capacities per volume (heat capacitances), c f , and c s , and thefluid velocity, v (Table 1). Table 1.
Material properties and constants used in calculations
Name Symbol Value or range Source
Fluid, solid density ρ f , ρ s (Lesher & Spera, 2015)Fluid specific heat capacity c p fluid c p solid c f . × J/(m K) c p fluid × ρ f Solid heat capacity per volume c s . × J/(m K) c p solid × ρ f Heat transfer coefficient k − to 10 W/m K this work (section 2.1)Fluid volume fraction φ v z Nu . . β A d − to 10 m (LeRoux et al., 2008) To non-dimensionalize the equations, we define the normalized relative tempera-ture of the fluid and solid, T (cid:48) f = ( T f − T ) / ∆ T and T (cid:48) s = ( T s − T ) / ∆ T , where T isreference temperature and ∆ T is a temperature perturbation (described below). We alsointroduce the dimensionless position, x (cid:48) = xk f /v , a dimensionless time, t (cid:48) = k s t , andthe heat capacitance ratio z = k s /k f = φc f / (1 − φ ) c s . The non-dimensional versionsof equations (1) and (2) are now (3) and (4): z ∂T (cid:48) f ∂t (cid:48) + ∂T (cid:48) f ∂x (cid:48) = − ( T (cid:48) f − T (cid:48) s ) (3) ∂T (cid:48) s ∂t (cid:48) = ( T (cid:48) f − T (cid:48) s ) (4)(see also Spiga and Spiga (1981)). It is clear that for a given temperature difference, ( T (cid:48) f − T (cid:48) s ), the behavior of this set of non-dimensional equations is governed by z (1 /z is thedimensionless fluid velocity). Analytic solutions for this set of equations have been de-rived for a number of limiting cases, particularly for large k (Spiga & Spiga, 1981; Kuznetsov,1994, 1995b, 1995a, 1996), and were used to benchmark the numerical calculations be-low. –4–anuscript submitted to Before we can interpret the results below, we consider the meaning of the assumedheat transfer coefficient, k , which is related to the individual solid and fluid heat trans-fer coefficients, k s = k/ ( c s (1 − φ )) and k f = k/ ( c f φ ). Since k f and k s have dimen-sions of inverse time, k represents the amount of heat transferred per unit time from fluidto solid, per unit volume, per unit difference in temperature (Schumann, 1929). Althoughthe non-dimensional system above (Eqns 3 and 4) is independent of k , the actual phys-ical behavior depends on the characteristic length ( v/k f ) and time (1 /k s ) scales, whichdepend on k .The factors that determine k are explored in Appendix A, but here we summarize:for a given fluid volume fraction, φ , k is strongly controlled by the length scale of chan-nelization, parameterized by the channel spacing d . To get an idea of the scale of chan-nelization at the LAB, we turn to structural, petrologic, and geochemical data from theLherz Massif suggesting that melt-rock interaction has driven refertilization of a harzbur-gite body into lherzolite (LeRoux et al., 2007, 2008). In the field, the lherzolite bodiesare separated from each other by distances of several tens of meters and this is also thespatial scale of isotopic disequilibrium between metasomatizing fluids and the harzbur-gite parent material (LeRoux et al., 2008). We take this as a proxy for the spatial sep-aration of fluid-rich channels and choose a broad range for the relevant spatial scale ofchannelization, d = 10 − to 10 m (10 cm to 100 m channel spacings). The correspond-ing range of the heat transfer coefficient (Table 1) is therefore k ≈ − to 10 W m − K − ,used in the models below. In the following, we shall consider modeling scenarios wherematerial properties are fixed and the flow regime is considered to be in the channelizedregime. For the first set of calculations, we assume the domain is initially at steady-statewith the fluid and solid in equilibrium, T s = T f = T . At t = 0, the temperature ofthe fluid entering at the inflow, x = 0, is perturbed so that, T f ( x = 0 , t ≥
0) = T +∆ T , introducing ∆ T as a temperature scale into the problem. This disturbs the initialsteady state and starting at t = 0 + , the fluid at the inflow is no longer in thermal equi-librium with material in the domain. Equations (3) and (4) may be solved for the ther-mal evolution of the fluid and solid subject to the initial and boundary condition: T (cid:48) s = T (cid:48) f = 0 initially and T (cid:48) f ( x = 0 , t ≥
0) = 1. The behavior of the nondimensional sys-tem is controlled by z , the dimensionless fluid velocity, which is a function of both ma-terial properties c s and c f and the fluid volume fraction, φ . Since the material proper-ties c s and c f are held constant, there is a unique mapping between z and φ , the fluidvolume fraction (we use φ = 0.01 to 0 .
20, corresponding to values of z = 0 . . k . The location of maximum disequilibrium (maximum T fluid − T solid ,and therefore the greatest heat exchange) lags behind the fluid front and progresses in-ward into the domain at a rate determined by the fluid velocity, fluid volume fraction,and material heat capacitances (Eqn 5; Appendix B). Models using a range of φ and v fluid values (Table 1) show that the disequilibrium front migrates at a rate given by V diseqm ≈ v fluid (cid:18) c f c s (cid:19) (cid:18) c f φc f φ + (1 − φ ) c s (cid:19) (5)These results are an extension of those in Kuznetsov (1994), where an analytic expres-sion for the migration rate is presented, in the limit that the degree of disequilibrium is –5–anuscript submitted to small. Although the rate of migration of the disequilibrium front is not a function of k ,the heat transfer coefficient, the characteristic width of this zone and the degree of dis-equilibrium within it are strong functions of k (Appendix B). Here we consider sinusoidal thermal pulses, introducing a timescale into the prob-lem, the period, τ , or the non-dimensional period τ k s = τ k/ (1 − φ ) c s , where 1 /k s isthe longest response timescale in the problem, associated with the thermal response ofthe solid. This class of models assume that the fluid entering the domain is hotter thanthe ambient initial temperature of the solid, but the thermal contrast varies sinusoidally.The time-varying inlet fluid temperature represents pulses of high temperature fluids ormelt, and the non-dimensional period controls the length scale, δ , over which thermaloscillations penetrate into the domain. For the material parameters in Table 1, and chan-nel spacing of d = 10 to 100 m leads to a response timescale 1 /k s ≈ /k s penetrate farther into the domain than short period oscillations(Figure 2). Therefore, periodic thermal perturbations that might represent melt infil-tration pulses lasting 10 to 10 yrs will be characterized by a region of sinusoidally vary-ing temperatures: a thermal re-working zone (TRZ) (blue curves in Figure 2a & b). Thewavelength of these oscillations is set by the period τ , λ = v fluid τ . The penetrationdistance of the thermal oscillations, δ , is the maximum width of the TRZ. At short times,when t (cid:28) δ/V diseqm , there is one zone of disequilibrium (the TRZ, bounded by the dis-equilibrium front). At longer times, the TRZ widens to a maximum width, δ , at time δ/V diseqm . When t (cid:29) δ/V diseqm t , there are two zones of disequilibrium, one station-ary at the inlet (the TRZ) and the migrating zone discussed above that moves at V diseqm (Eqn 5; red curves in Figure 2a & b).To illustrate the oscillatory nature of temperatures inside the TRZ, we consider tem-perature vs. time paths within the domain at varying distances from the inlet (Figure2c). The amplitude of the oscillations decay with distance, but at each point in the TRZ,once the oscillations are established, the amplitude is constant in time (Figure 2c). Aswe might expect, the maximum width of the TRZ, δ is set both by the non-dimensionaloscillation period, τ /t s and by the heat transfer coefficient (Spiga & Spiga, 1981), δ = (cid:18) c f φv fluid k (cid:19) (cid:18) ( τ /t s ) π (cid:19) (6)noting that k ∼ d − (Appendix A), and t s = 1 /k s = c s (1 − φ ) /k , the expressionabove suggests that, for fixed v fluid , δ ∼ ( τ /d ) as confirmed by the numerical results(Figure 2d). The models presented here are highly idealized and therefore limited in their rep-resentation of the complexities of deformation and fluid-rock interactions within the Earth.In particular, we abstract the effective thermal properties and the geometry of the fluid-solid interface into a single number, the heat transfer coefficient, k , which is strongly con-trolled by the channel spacing, d . Sinuosity and other aspects of the geometry of chan-nelization are abstracted and the details of processes at and below the scale of an av-erage channel spacing, d , are ignored. Instead, the focus here is on the effective behav-ior at mesoscopic spatial scales (cid:29) d . Even at these scales, we ignore spatial variationsin transport, including variability in the fluid volume fraction, fluid velocity, and effec-tive heat transfer coefficient. Additionally we have ignored time-dependent variabilityin transport, e.g., feedbacks due to possible phase changes during disequilibrium heat-ing/cooling, which would affect the geometry of the fluid-solid interface (re. Keller and –6–anuscript submitted to (a)(b) δ δ δ = width, thermal reworking zone (a) (c) (d) Figure 2. (a) Normalized temperature profiles, T s / ∆ T (dashed) and T f / ∆ T (solid) for acalculation with fluid velocity v = 1 m/yr, at time t = 1 Kyr. Results are shown for two caseswith the same fluid velocity, φ , and channel spacing d indicated, but with different normalizedperiods, τ k s = 50 (red), 150 (blue), for the thermal perturbation. For the chosen parameters, t s = 1 yr. The amplitude of spatial oscillations decreases over a decay scale δ , as indicated. (b) Thedegree of disequilibrium ( T (cid:48) f − T (cid:48) s ) is also oscillatory over a thermal reworking zone of width δ .(c) Temperature-time paths within x < δ plotted at different distances from the inlet for the casewhere τ /t s = 150 in (a) and (b), illustrating temporal oscillations at each location within thethermal reworking zone. (d) Width of the thermal reworking zone, δ , as a function of oscillationperiod τ and channel spacing d as indicated. Dashed lines (slope 2) are the expected analyticscaling in Eqn 6 and the squares indicate numerically derived values of δ obtained by fitting anexponential decay to the envelope of the ( T (cid:48) f − T (cid:48) s ) oscillations in the TRZ, e.g., in (b). Thinhorizontal lines are at δ = 1 and 10 km. Suckale). Finally, the models presented here are 1D and ignore the 3D nature of rela-tive motion between fluid and solid even on the mesoscale ( (cid:29) d ).Given these limitations, the 1D models are presented here as a useful way to framefirst-order questions and develop arguments related to the consequences of disequilib-rium thermal transport dominated by downstream effects in the direction of fluid flow.First, lets assume that the model domain is analogous to the lowermost lithosphere, wherewe expect melt or fluid transport to be channelized (e.g., Schmeling et al., 2018, Figure –7–anuscript submitted to x < V fluid = 0.1 to 1 m/yr max. width, δ region of diseqm heating d TRZ width
10 km m V diseqm =10 km/Myr Initial LABTop of TRZ at t=0.1 to 1 Myr V fluid = 0.1 - 1 m/yr Melt pulse timescale τ =10,000 yr Figure 3.
Cartoon illustrating implications for a thermal re-working zone (TRZ) that formsa modified layer at the lowermost CLM as a result of disequilibrium heating. For episodic melt-infiltration into channels of spacing d = 100 m, with period ∼
10 Kyr, the TRZ is characterizedby upward-decreasing degree of disequilibrium (indicated by the color) and is δ ∼
10 km wideafter 0.1 to 1 Myr, for v fluid = 0 . gion (e.g. a decompaction layer, Holtzman & Kendall, 2010), whereas the domain x > x = 0),melt-infiltration into the lithosphere may be episodic, controlled by timescales associ-ated with transport (e.g., Scott & Stevenson, 1984; Wiggins & Spiegelman, 1995), pro-cesses of fracturing and crystallization in a diking 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, es-timated using the heat transfer coefficient, may be a significant portion of the heat bud-get at the LAB and the lowermost CLM, (2) a material-dependent velocity scale asso-ciated with transient disequilibrium heating, and (3) the existence of a thermal rework-ing zone (TRZ) associated with spatio-temporally varying disequilibrium heat exchange.We discuss each of these within the context of episodic melt-infiltration into the base ofthe CLM in an intra-plate setting. Specifically, an intra-plate tectonic setting such asBasin and Range province of the western US where deformation and 3D melt-rock in-teraction may be simplified so that we consider the effects of dominantly vertical heattransport within a slowly deforming lithosphere. i. Disequilibrium heating and the heat budget at the LAB. The relativeimportance of disequilibrium heating at the LAB may be established by considering theeffective heat transfer coefficient, k , and the factor which most strongly controls it, namelythe average spacing of channels, d . For the material parameters in Table 1, and chan-nel spacing of d = 1 to 100 m, k is in the range k ≈ − to 10 W m − K − (AppendixA). Physically, k corresponds to fluid-solid heat transfer per unit time, per unit volume,per unit difference in temperature (Schumann, 1929). Therefore, for a 100 K excess tem-perature of the infiltrating melt, disequilibrium heating might contribute around 10 − to 10 W m − to the heat budget at the LAB. This is a conservative estimate, given that –8–anuscript submitted to the temperature difference between magma and the surrounding material may be any-where from 100 to 1000 K (e.g., in crust; (Lesher & Spera, 2015)). (While an LAB ex-cess temperature of the melt of about 100 K is reasonable, but note that plume excesstemperatures 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 heatbudget due to crystallization of melt in channels may be estimated using scaling argu-ments made in Havlin et al. (2013). Assuming that melt and rock are in equilibrium, Havlinet al. (2013) estimate that the heat released by a crystallization front would contributearound ρHS dike , where ρ is the melt density, H is the latent heat of crystallization, S dike is a volumetric flow rate out of a decompacting melt-rich LAB boundary layer due todiking. For a representative porosity of φ = 0 . S dike ≈ × − m /s. Taking ρ = 3000 kg m , and H = 3 × J/kg, theheat source due to the moving crystallization front would be around 10 W for each 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 power per unit volume generateddue to crystallization. For example, assuming dike heights of about 10 m and dike spac-ing large enough for non-interacting dikes (as estimated by Havlin et al. (2013), a poros-ity of 0.1 would require a dike spacing of ∼ m), the heat source due to a crystal-lizing dike boundary layer would be < − W/m (per unit length along strike). Thesearguments corroborate the idea that disequilibrium heating during melt-rock interactioncould be a significant portion of the heat budget at the LAB as compared to other ex-pected processes, such as heating due to crystallization of melt in channels. ii. Progression of a disequilibrium heating zone/front at a rate V diseqm . The importance of a reduced velocity for the migration rate of either a disequilibriumfront (Figure B1a,b,c) or widening of a thermal reworking zone (Figure 2b) is that thislimits the rate at which the lowermost CLM may be modified by thermal disequilibrium.It is important to note that this reduced velocity (Eqn 5) is independent of temperaturecontrast between the CLM and infiltrating melt and depends only on the relative fluidvolume fraction, fluid velocity and material properties. Assuming a fluid volume frac-tion of 1 to 10 % at the LAB, and material properties in Table 1, we would expect that V diseqm would be around 1 to 10 % of the fluid velocity (see Appendix B and FiguresB1, 3). For relative fluid velocity in the range of 0.01 to 1 m/yr, we would predict thatdisequilibrium heating front at the LAB would migrate upward at a rate of ≈ km/Myr, which is comparable to rates of CLM thinning predicted by heating due to theupward motion of a dike boundary layer (1 to 6 km/Myr in Havlin et al. (2013)). Inter-estingly, an upward-moving disequilibrium heating zone with V diseqm ≈ km/Myrbrackets the rate of upward migration of the LAB at 10-20 km/Myr inferred from thepressure and temperature of last equilibration of Cenozoic basalts in the Big Pine vol-canic Field in the western US (Plank & Forsyth, 2016). An implication of the modelshere, therefore, is that disequilibrium heating may produce lithosphere modification atgeologically-relevant spatial and temporal scales provided that the melt velocity in chan-nels at the LAB is on the order of 10 − to 10 − m/yr (Figure 3). iii. Thermal reworking zone (TRZ). A key result that may be relevant to theevolution of the LAB is that episodic infiltration of melts that are hotter than the sur-rounding CLM would lead to a long-lived region of disequilibrium heating within a ther-mal reworking zone or TRZ. Although the models suggest that the TRZ would undergoa phase of transient widening (at a rate given by Eqn 5), it will reach a steady-state width δ that should scale as δ ∼ (cid:2) φv fluid d τ (cid:3) where d is a characteristic scale of channel-ization and τ is a timescale associated with the episodicity of melt-infiltration (Figures2d and 3). This scaling gives us a way to conceptualize the modification of the lower-most CLM as a zone that encompasses a variable thickness TRZ, depending on v fluid and on the timescale of melt-infiltration (Figure 3). Regions where the timescale of episodicmelt-infiltration is longer are predicted to have a thicker zone of modification at the LAB. –9–anuscript submitted to For example, for a channel spacing of d = 10 m, disequilibrium heating by repeatedmelt pulses that last around 10 Kyrs implies a maximum thickness of roughly 10 km forthe zone of modification (Figure 2d). In this scenario, the TRZ grows to this maximumwidth over a timescale governed by δ/V diseqm ; for V diseqm = 10 km/Myr, which cor-responds to melt velocity of roughly 0.1 m/yr (see (ii) above), the 10 km wide TRZ wouldbe established within about 1 Myr (Figure 3), comparable to that icomparable to ratesof CLM modification inferred from observations in Plank and Forsyth (2016). In summary, the spatial and temporal scales associated with the establishment ofthe TRZ are comparable to those for CLM modification inferred from geochemical andpetrologic observations intra-plate settings, e.g., the western US (Plank & Forsyth, 2016).These scaling arguments lead to the idea that perhaps the TRZ represents a zone of ther-mal 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 hotterthan the ambient material. The dynamic evolution of the LAB during episodic pulsesof melt-infiltration is beyond the scope of the simple models above (which assume a sta-tionary, undeforming matrix). However, assuming mantle material that obeys a temper-ature and pressure-dependent viscosity scaling relation such as in Hirth and Kohlstedt(2003), at an LAB depth of about 75 km where we assume that the ambient mantle iscooler than the dry solidus (e.g., 1100 o C + 3.5 o C/km; Plank and Forsyth (2016)), wewould expect a significant viscosity reduction during heating (e.g., factor ≈ /
62 for atemperature increase of 100 K). This effect is weaker, but still important for a deeperLAB; e.g, at 125 km depth, the viscosity reduction would be a factor ≈ /
18 for a tem-perature increase of 100 K.The models above suggest that disequilibrium heating may contribute more than10 − W/m to the heat-budget at the LAB and, for fluid velocity of 0.1 to 1 m/yr in chan-nels that are roughly 10 m apart, a 10 km wide TRZ may be established within 1 Myr.Disequilibirum heating, therefore, is associated with a timescale that may be relevantfor regional tectonics in intra-plate settings. Interestingly, geochemical evidence from Ceno-zoic basalts from the western US, particularly space-time variations in volcanic rock Ta/Thand Nd isotopic compositions suggest that the timescale of modification and removal ofthe lowermost CLM is on the order of 10 Myrs (Farmer et al., 2020). These authors ar-gue that the observed transition from low to intermediate to high Ta/Th ratios indicatesa change from: arc/subduction-related magmatism, to magmatism associated with in-situ melting of a metasomatized CLM (the “ignimbrite flare-up”), to magmatism dueto decompression and upwelling after removal of the lowermost CLM. If correct, theseinterpretations and observations are promising and provide an important avenue for ex-ploring the role of thermal and chemical disequilibrium during melt-rock interaction anddestabilization of the CLM in an intra-plate setting. Disequilibrium heating during melt-infiltration may be an important process for modifying the lowermost CLM and may playa role in the rheologic weakening that may be needed to mobilize and possibly removethe lowermost CLM.
Appendix A Heat transfer coefficient, k The factors that determine k can be illustrated by considering that the fluid-solidheat transfer rate must depend on the geometry of the inter-phase interface and also onthe effective thermal conductivity of the porous medium. Although the geometry of thefluid-solid interface may be complex, this model considers one aspect of it: the specificfluid-solid interface surface area (contact area per unit volume), a sf , which is a functionof the length-scale of channelization in the solid. In the porous flow case for example,if the solid matrix is made of spheres with an average particle diameter d , then the spe- –10–anuscript submitted to cific area for a grain is S = 6 /d , so a sf = S (1 − φ ) = 6(1 − φ ) /d (Dullien, 1979).This sets a limit for channels, where we shall assume that the specific surface area is a sf ∼ A (1 − φ ) /d , where A is a number that is between 2 (as for a single cylindrical channelwith small volume fraction φ ) and 6 (Dixon & Cresswell, 1979). Whereas the specific con-tact area is a geometric factor, the effective conductivity of the medium depends on theNusselt number, N u . Theoretical arguments in Dixon and Cresswell (1979) show thatthe effective thermal conductivity may be written in terms of the individual fluid andsolid thermal conductivities λ f and λ s (basically taking the fluid and solid in parallel):1 C eff = (cid:20) N uλ f − βλ s (cid:21) (A1)where β = 10 for spherical matrix grains, 8 for cylinders, and 6 for slabs (Dixon & Cress-well, 1979). For the rest of this work, we will take the range of β to be 6 to 10 represent-ing the highly-channelized vs porous flow end-member geometries. For slow flows (Reynoldsnumber Re (cid:28) N u ranges from 0 . . k is an ef-fective “conductance” C eff /d , so that k = C eff a sf /d , k = 1 d (cid:20) N uλ f − βλ s (cid:21) − A (1 − φ ) d (A2)a product of a material-dependent quantity and a geometry-dependent quantity. Turn-ing now to physical properties relevant to the transport of melts through the lithosphere,using a reasonable conductivity for basaltic magma of λ f = 1 W/(m K) (Lesher & Spera,2015), λ s = 2 . N u = 0 . . β = 6, 8, or10, and A = 6, we find that the effective conductivity C eff , a sf and k are within theranges shown in Figure A1. (Note that choosing A ≈
2, for a more channelized geom-etry, would change k by less than an order of magnitude.) Specifically, Eqn (A2) showsthat k ∼ d − and strongly decreases with increasing spatial scale of channels charac-terizing the fluid-solid interaction; Figure A1c. Appendix B Response to a step-function
As fluid with a perturbed temperature enters at x = 0, the solid heats up whilethe fluid cools (for positive perturbation ∆ T ). The perturbation “front”, the farthestextent of the fluid that entered x = 0 at t = 0 with perturbed temperature, is alwaysat x front = vt (or x (cid:48) front = t (cid:48) /z ; where the dimensionless velocity of the fluid is 1 /z ).The response to a step-function is essentially a transient disequilibirum front, travelinginward at V diseqm , behind which the solid and fluid equilibrate at the new inlet temper-ature.As we might expect, the behavior of the nondimensional system (Eqns 3 and 4) isgoverned by z , the heat capacitance ratio. In response to a step-function increase in thefluid temperature at the inlet, temperature profiles within the domain exhibit a tran-sition or disequilibrium zone, lagging behind the fluid front (Figure B1). Ahead of thisdisequilibrium zone, the fluid and solid are in equilibrium at the initial ambient temper-ature of the solid, T (cid:48) s = T (cid:48) f = 0. Behind this zone, the fluid and solid are in equilib-rium at the incoming fluid temperature, T (cid:48) s = T (cid:48) f = 1.xFollowing an initial lag time (when the maximum disequilibrium is at x (cid:48) = 0), thedisequilibrium zone migrates inward migration at a steady speed, a fixed fraction of thefluid velocity, v (Figure B1c and B1d). The ratio of the migration rate of the disequi-librium zone to the fluid velocity is controlled by the heat capacity ratio, z , and there-fore by the fluid volume fraction, φ (Figure B1). In the limit that the fluid and solid arenearly in equilibrium ( T (cid:48) f − T (cid:48) s ≈ –11–anuscript submitted to (a) (b)(c) Figure A1. (a) Effective thermal conductivity, C eff , as a function of Nusselt number, (b)geometric factor, a sf , as a function of transport length scale in the solid, d , and (c) fluid-solidheat transfer coefficient, k , as a function of transport length scale in the solid, d . For a fixed d value, 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. on √ t (cid:48) and the zone of disequilibrium migrates at speed vc f / ( φc f + (1 − φ ) c s ). Our mod-els show that, when there is significant disequilibrium, the zone of disequilibrium migrateswith a rate given by Eqn 5 (Figure B1d), which does not depend on k , the heat trans-fer coefficient.We illustrate the dependence on k for the specific case where we take the fluid ve-locity v =1 m/yr, and consider channel spacings, d = 50 to 150 m, which correspondto k = 5 × − and 5 × − W m − K − , respectively (Figure B2).The migration of the disequilibrium front may be thought of as the motion of thelocus of maximum heating, which moves at speed V diseqm ≈ v/
10, for φ = 0 . D , de-cays as 1 / √ t (e.g., Kuznetsov, 1994) and is controlled by the heat transfer coefficient, k . D scales as 1 / √ k and therefore is a linear function of d , the channel spacing (FigureB2b). This transient is apparent in temperature-time paths (Figure B2c and B2d), wherethe approach to steady-state occurs on a timescale governed by φv fluid /k . Acknowledgments
This work benefited from conversations with B. Holtzman, L. Farmer, and R. Carlsonregarding western US evolution, and B. Oliveira on melt-rock interaction. The codes (writ-ten in Matlab) used for the models presented are available upon request. –12–anuscript submitted to
Figure B1. (a) Normalized temperature as a within the fluid (solid lines) and solid (dashedlines) at different values of the dimensionless time, t (cid:48) , as indicated. The colors represent two dif-ferent fluid volume fractions, φ , and therefore different z . (b) Normalized temperature differencebetween fluid and solid as a function of dimensionless position, shown for the cases consideredin (a). (c) The same profiles as in (b), but now plotted as a function of position normalized bythe fluid front location. The fluid front is always at x (cid:48) front = t (cid:48) /z ; stationarity of the disequilib-rium zone in this plot indicates that the disequilibrium zone migrates at a constant, z-dependentfraction of the fluid velocity. (d) Normalized migration rate of the zone of disequilibrium as afunction of fluid volume fraction, φ . Red dot is for φ ≈ References
Bo, T., Katz, R. F., Shorttle, O., & Rudge, J. F. (2018). The melting column as afilter of mantle trace-element heterogeneity.
Geochemistry, Geophysics, Geosys-tems , , 4694–4721.Carlson, R. W., Irving, A. J., Schulzec, D. J., & Jr, B. C. H. (2004). Timing ofprecambrian melt depletion and phanerozoic refertilization events in the litho-spheric mantle of the wyoming craton and adjacent central plains orogen. Lithos , (453-472).Cashman, K. V., Sparks, R. S. J., & Blundy, J. D. (2017). Vertically extensive andunstable magmatic systems: A unified view of igneous processes. Science , ,1280.Dixon, A. G., & Cresswell, D. L. (1979). Theoretical prediction of effective heattransfer parameters in packed beds. AIChE Journal , (4).Dullien, F. A. L. (1979). Porous media: Fluid transport and pore structure . Aca-demic Press (New York).Farmer, G. L., Fritz, D., & Glazner, A. F. (2020). Identifying metasomatized con-tinental lithospheric mantle involvement in cenozoic magmatism from ta/thvalues, southwestern north america.
Geochemistry, Geophysics, Geosystems , , 10.1029/2019GC008499.Gao, S., Rudnick, R. L., Carlson, R. W., McDonough, W. F., & Liu, Y.-S. (2002).Re-os evidence for replacement of ancient mantle lithospherebeneath the north –13–anuscript submitted to (a)(b) wD Figure B2. (a) Normalized temperature profiles, T s / ∆ T (dashed) and T f / ∆ T (solid) for acalculation with fluid velocity v = 1 m/yr, at times t = 1, 12 .
5, and 25 Kyr. The tempera-ture profiles transition between the incoming fluid temperature (=1, on the left) and the initialambient temperature (=0, on the right). Results are shown for two cases with different channelspacings, d , corresponding to different heat transfer coefficients, k , as indicated. The transitionregion (e.g., highlighted in gray at t =12 . w , that is larger forsmaller k (large d ) and increases over time. (b) The degree of disequilibrium within the transitionzone is characterized by the maximum difference, D , between the solid and fluid temperatureprofiles. D is greater for smaller k and decreases as a function of time. (c) and (d), Normalized(c) solid and (d) fluid temperature as a function of dimensionless time since first contact with theperturbation front for the results in blah. china craton. Earth and Planetary Science Letters , , 307-322.Handley, D., & Heggs, P. J. (1968). Momentum and heat transfer mechanisms inregular shaped packings. Trans. Inst. Chem. Engrs. , , T251.Havlin, C., Parmentier, E., & Hirth, G. (2013). Dike propagation driven by meltaccumulation at the lithosphere-asthenosphere boundary. Earth and PlanetaryScience Letters , , 20-28. doi: 10.1016/j.epsl.2013.06.010Hirth, G., & Kohlstedt, D. L. (2003, Dec). Rheology of the upper mantle and themantle wedge: A view from the experimentalists. Geophysical Monograph Se-ries , , 83–105.Holtzman, B. K., & Kendall, J.-M. (2010, Dec). Organized melt, seismic anisotropy, –14–anuscript submitted to and plate boundary lubrication. Geochem. Geophys. Geosyst. , , Q0AB06.doi: 10.1029/2010GC003296Hopper, E., Gaherty, J. B., Shillington, D. J., Accardo, N. J., Nyblade, A. A., Holtz-man, B. K., . . . Mbogoni, G. (2020). Preferential localized thinning of litho-spheric mantle in the melt-poor malawi rift. Nature Geoscience , , 584-589.Keller, T., & Suckale, J. (2019). A continuum model of multi-phase reactive trans-port in igneous systems. Geophys. J. Int. , , 185-222.Kuznetsov, A. V. (1994). An investigation of a wave of temperature difference be-tween solid and fluid phases in a porous packed bed. International Journal ofHeat and Mass Transfer , (18), 3030-3033.Kuznetsov, A. V. (1995a). An analytical solution for heating a two-dimensionalporous-packed bed by a non-thermal equilibrium fluid flow. Applied ScientificResearch , (1), 83-93.Kuznetsov, A. V. (1995b). Comparison of the waves of temperature difference be-tween the solid and fluid phases in a porous slb and in a semi-infinite porousbody. International Communications in Heat and Mass Transfer , (4),499-506.Kuznetsov, A. V. (1996). Analysis of heating a three-dimensional porous bed utiliz-ing the two energy equation model. Heat and Mass Transfer , (3), 173-178.Lamb, S., Moore, J. D. P., Smith, E., & Stern, T. (2017). Episodic kinematics incontinental rifts modulated by changes in mantle melt fraction. Nature , .LeRoux, V., Bodinier, J., Tommasi, A., Alard, O., Dautria, J., Vauchez, A., &Riches, A. (2007). The lherz spinel lherzolite: Refertilized rather than pristinemantle. Earth and Planetary Science Letters , , 599–612.LeRoux, V., Bodinier, J.-L., Alard, O., & S.Y. O’Reilly, W. G. (2008). Isotopic de-coupling during porous melt flow: A case-study in the lherz peridotite. Earthand Planetary Science Letters , , 76-85.Lesher, C. E., & Spera, F. J. (2015). Chapter 5 - thermodynamic and transportproperties of silicate melts and magma. In (p. 113-141). Elsevier.Menzies, M., Xu, Y., Zhang, H., & Fan, W. (2007). Integration of geologygeophysicsand geochemistry: A key to understanding the north china craton. Lithos , ,1–21.Oliveira, B., Afonso, J., Zlotnik, S., & Diez, P. (2018). Numerical modelling of mul-tiphase multicomponent reactive transport in the earth’s interior. Geophys. J.Int. , (1), 345-388.O’Reilly, S. Y., Griffin, W. L., Djomani, Y. H. M., & Paul. (2001). Are lithospheresforever? Tracking changes in the subcontinental lithospheric mantle throughtime. GSA Today , (4), 4–10.Plank, T., & Forsyth, D. W. (2016). Thermal structure and melting conditions inthe mantle beneath the basin and range province from seismology and petrol-ogy. Geochemistry Geophysics Geosystems , (1312-1338).Roy, M., Gold, S., Johnson, A., Orozco, R. O., Holtzman, B. K., & Gaherty, J.(2016). Macroscopic coupling of deformation and melt migration at continentalinteriors, with applications to the colorado plateau. Journal of GeophysicalResearch , , doi:10.1002/2015JB012149.Schmeling, H., Marquart, G., & Grebe, M. (2018). A porous flow approach to modelthermal non-equilibrium applicable to melt migration. Geophysical Journal In-ternational , (119-138).Schumann, T. E. (1929). Heat transfer: a liquid flowing through a porous prism. Journal of the Franklin Institute , (3), 405–416.Scott, D. R., & Stevenson, D. J. (1984). Magma solitons. Geophys Res Lett , (11),1161-1164.Spiga, G., & Spiga, M. (1981). A rigorous solution to a heat transfer two-phasemodel in porous media and packed beds. International Journal of Heat andMass Transfer , , 355-364. –15–anuscript submitted to Wallner, H., & Schmeling, H. (2016). Numerical models of mantle lithosphereweakening, erosion and delamination induced by melt extraction and emplace-ment.
International Journal of Earth Science (Geologische Rundschau) , ,1741-1760.Wang, H., van Hunen, J., & Pearson, D. G. (2015). The thinning of subcontinentallithosphere: The roles of plume impact and metasomatic weakening. Geochem-istry Geophysics Geosystems , , doi:10.1002/ 2015GC005784.Wenker, S., & Beaumont, C. (2017). Can metasomatic weakening result in the rift-ing of cratons? Tectonophysics , 1940-1951.Wiggins, C., & Spiegelman, M. (1995). Magma migration and magmatic solitarywaves in 3-d.
Geophys Res Lett , (10), 1289-1292.(10), 1289-1292.