Articles | Volume 13, issue 6
Research article
27 Jun 2022
Research article |  | 27 Jun 2022

Thermal non-equilibrium of porous flow in a resting matrix applicable to melt migration: a parametric study

Laure Chevalier and Harro Schmeling

Fluid flow through rock occurs in many geological settings on different scales, at different temperature conditions and with different flow velocities. Depending on these conditions the fluid will be in local thermal equilibrium with the host rock or not. To explore the physical parameters controlling thermal non-equilibrium, the coupled heat equations for fluid and solid phases are formulated for a fluid migrating through a resting porous solid by porous flow. By non-dimensionalizing the equations, two non-dimensional numbers can be identified controlling thermal non-equilibrium: the Péclet number Pe describing the fluid velocity and the porosity ϕ. The equations are solved numerically for the fluid and solid temperature evolution for a simple 1D model setup with constant flow velocity. This setup defines a third non-dimensional number, the initial thermal gradient G, which is the reciprocal of the non-dimensional model height H. Three stages are observed: a transient stage followed by a stage with maximum non-equilibrium fluid-to-solid temperature difference, ΔTmax, and a stage approaching the steady state. A simplified time-independent ordinary differential equation for depth-dependent (TfTs) is derived and solved analytically. From these solutions simple scaling laws of the form Tf-Ts=fPe,G,z are derived. Due to scaling they do not depend explicitly on ϕ anymore. The solutions for ΔTmax and the scaling laws are in good agreement with the numerical solutions. The parameter space PeG is systematically explored. Three regimes can be identified: (1) at high Pe (>1/G) strong thermal non-equilibrium develops independently of Pe, (2) at low Pe (<1/G) non-equilibrium decreases proportional to decreasing PeG, and (3) at low Pe (<1) and G of the order of 1 the scaling law is ΔTmaxPe. The scaling laws are also given in dimensional form. The dimensional ΔTmax depends on the initial temperature gradient, the flow velocity, the melt fraction, the interfacial boundary layer thickness, and the interfacial area density. The time scales for reaching thermal non-equilibrium scale with the advective timescale in the high-Pe regime and with the interfacial diffusion time in the other two low-Pe regimes. Applying the results to natural magmatic systems such as mid-ocean ridges can be done by estimating appropriate orders of Pe and G. Plotting such typical ranges in the PeG regime diagram reveals that (a) interstitial melt flow is in thermal equilibrium, (b) melt channeling such as revealed by dunite channels may reach moderate thermal non-equilibrium with fluid-to-solid temperature differences of up to several tens of kelvin, and (c) the dike regime is at full thermal non-equilibrium.

1 Introduction

Fluid flow through rock occurs in many geological settings on different scales, at different temperature conditions and with different flow velocities. Depending on these conditions the fluid will be in local thermal equilibrium with the host rock or not. On a small scale, e.g., grain scale, usually thermal equilibrium is valid. Examples include melt migration through a porous matrix in the asthenosphere or in crustal magmatic systems at super-solidus temperatures (e.g., McKenzie, 1984), groundwater or geothermal flows in sediments or cracked rocks (e.g., Verruijt, 1982; Furbish, 1997; Woods, 2015), or hydrothermal convection in the oceanic crust (e.g., Davis et al., 1999; Harris and Chapman, 2004; Becker and Davies, 2004). On a somewhat larger scale local thermal equilibrium may not always be reached. Examples of such flows include melt migration in the mantle or crust at temperatures close to or slightly below the solidus where melt may be focused and migrates through systems of veins or channels (Kelemen et al., 1995; Spiegelman et al., 2001). Within the upper oceanic crust water may also migrate through systems of vents or channels (Wilcock and Fisher, 2004). At even larger scales and at sub-solidus conditions, magma rapidly flows through propagating dikes or volcanic conduits (e.g., Lister and Kerr, 1991; Rubin, 1995; Rivalta et al., 2015) and is locally at non-equilibrium with the host rock.

Heat transport associated with most of such flow scenarios is usually described by assuming thermal equilibrium between the fluid and solid under slow flow conditions (e.g., McKenzie, 1984). Alternatively, for more rapid flows such as melts moving in dikes through a cold elastic or visco-elasto-plastic ambient rock, the fluids are assumed to be isothermal (e.g., Maccaferri et al., 2011; Keller et al., 2013). However, on a local scale of channel or dike width, thermal interaction between rising hot magma and cold host rock is important and may lead to effects such as melting of the host rock and freezing of the magma with important consequences for dike propagation and the maximum ascent height (e.g., Bruce and Huppert, 1990; Lister and Kerr, 1991; Rubin, 1995). Clearly, in such rapid fluid flow scenarios melt is not in thermal equilibrium with the ambient rock.

Thus, there exists a transitional regime, which, for example, may be associated with melt focusing into pathways where flow is faster and thermal equilibrium might not be valid anymore. In such a scenario it might be possible that channelized flow of melt might penetrate deeply into sub-solidus ambient rock, and thermal non-equilibrium delays freezing of the ascending melts and promotes initiation of further dike-like pathways. Indeed, for mid-oceanic ridges compositional non-equilibrium has proven to be of great importance for understanding melt migration and transport evolution (Aharonov et al., 1995; Spiegelman et al., 2001). Thus, it appears plausible that in cases of sufficiently rapid fluid flow, e.g., due to channeling or fracturing, thermal non-equilibrium may also become important. Describing this non-equilibrium macroscopically, i.e., on a scale larger than the pores or channels, is the scope of this paper.

While the physics of thermal non-equilibrium in porous flow is well studied in more technical literature (e.g., Schumann, 1929; Spiga and Spiga, 1981; Kuznetsov, 1994; Amiri and Vafai, 1994; Minkowycz et al., 1999; Nield and Bejan, 2006; de Lemos, 2016), so far it has attracted only little attention in the geoscience literature, but see Schmeling et al. (2018) and Roy (2020). The basic approach in all these studies is the decomposition of the heat equation for porous flow into two equations, one for the solid and one for the migrating fluid. The key parameter for thermal non-equilibrium is a heat exchange term between fluid and solid, which appears as a sink in the equation for the fluid and as a source in the equation for the solid. Usually, this heat exchange term is assumed proportional to the local temperature difference between fluid and solid (Minkowycz et al., 1999; Amiri and Vafai, 1994; de Lemos, 2016; Roy, 2020). However, Schmeling et al. (2018) showed that in a more general formulation the heat exchange term depends on the complete thermal history of the moving fluid through the possibly also moving solid. Here we will follow the common assumption and use the local temperature difference formulation. While Schmeling et al. (2018) showed that the magnitude of thermal non-equilibrium essentially depends on the flow velocity, or more precisely on the Péclet number, here we will more generally explore the parameter space.

While thermal non-equilibrium of an arbitrary porous flow system depends on many parameters, our approach is to reduce the complexity of the system and systematically explore the non-dimensional parameter space. It will be shown that only two non-dimensional parameters control thermal non-equilibrium in porous flow, namely the Péclet number and the porosity. In our simple 1D model setup with constant flow velocity a third non-dimensional number, the non-dimensional initial thermal gradient G, is identified, which is equal to the reciprocal non-dimensional model height H=1/G. The non-dimensionalization allows application of the results to arbitrary magmatic or other systems. The aim is to derive scaling laws that allow an easy determination of whether thermal equilibrium or non-equilibrium is to be expected and quantitatively to estimate the maximum temperature difference between fluid and matrix. The results will be applied, the magmatic system of a mid-ocean ridge setting.

2 Governing equations and model setup

2.1 Heat conservation equations

We start with considering a general two-phase matrix–fluid system with variable properties and solid and fluid velocities and subsequently apply simplifications. The two phases are incompressible, and we assume local thermal non-equilibrium conditions; i.e., the two phases exchange heat. The equations for conservation of energy of this system are given by de Lemos (2016), for example. Assuming constant pressure, the conservation of energy of the fluid phase is given by

(1) c p , f ϕ ρ f T f t + ϕ ρ f v f T f = ϕ λ f T f - Q fs .

For the definition of all quantities, see Table 1. Equation (1) can be rearranged into

(2) c p , f T f ϕ ρ f t + ϕ ρ f T f t + T f ϕ ρ f v f + ϕ ρ f v f T f = ϕ λ f T f - Q fs .

Mass conservation for the fluid phase is given by

(3) ρ f ϕ t + ρ f ϕ v f = 0 .

Inserting Eq. (3) into Eq. (2), conservation of energy for the fluid phase becomes

(4) c p , f ρ f ϕ T f t + v f T f = ϕ λ f T f - Q fs .

In a similar way, the conservation of energy of the solid phase is given by

(5) c p , s ρ s 1 - ϕ T s t + v s T s = 1 - ϕ λ s T s + Q fs ,

which, assuming that vs=0, is further simplified:

(6) c p , s ρ s 1 - ϕ T s t = 1 - ϕ λ s T s + Q fs .

The term Qfs in the fluid and solid heat conservation equations is the interfacial heat exchange term between the two phases (fluid and solid). In general, it depends on the local thermal history of the two phases and the history of the heat exchange (Schmeling et al., 2018). In a simplification it can be written as a combination of the interfacial area density S with the dimension [m−1], the interfacial boundary layer thickness δ, the effective thermal conductivity λeff, and the temperatures of the two phases:

(7) Q fs = S λ eff δ ( T f - T s ) .

In general, the term δ is time dependent. Schmeling et al. (2018), however, provide evidence that taking an appropriate constant value for δ (depending on fluid velocity) gives a good approximation of Qfs and allows for reasonable modeling of temperature evolution with time. In most of the following parametric study, we use this simplification for δ by assuming it is constant with time.

Table 1Symbols, their definition, and physical units used in this study.

Download Print Version | Download XLSX

2.2 Scaling and non-dimensionalization

Non-dimensionalization is useful for interpreting models involving a large number of parameters. It usually helps reduce the number of parameters and identifies non-dimensional parameters that control the evolution of the system. We write the two energy conservation equations in a non-dimensional form, using

(8) T = Δ T 0 T , t = t 0 t , v = v f 0 v , x , y , z = L x , y , z ,

where ΔT0 is the macroscopic scaling temperature difference of the system, i.e., the initial temperature difference between top and bottom; x,y, and z are distance; vf0 is the scaling fluid velocity; L is the scaling length

(9) L = ϕ 0 1 - ϕ 0 δ S ,

with ϕ0 as a scaling porosity; and t0 is the scaling time based on the diffusion time over the length L,

(10) t 0 = L 2 / κ 0

(see Table 1 for definitions). Primed quantities are non-dimensional.

Introducing the fluid-filled pore width df and the solid width ds, which may be the grain size or distance between fluid channels, the interfacial area density S scales with

(11) S = c ϕ 0 d f

for melt channels, tubes, pockets for all melt fractions, and melt films at small melt fractions, while S scales with

(12) S = c s 1 - ϕ 0 d s

for melt channels, films, and suspensions at all melt fractions. Here c is a geometrical constant of the order of 2 for melt channels, 4 for melt tubes, 6 for melt pockets, and 2 for melt films at small melt fractions. It should be emphasized that Eqs. (11) and (12) are different ways of calculating the same S. The geometrical constant cs is of the order of 2 for melt channels and 6 for melt films or suspensions. Thus, the scaling time and scaling length can also be written as

(10a) t 0 = 1 - ϕ 0 d f δ c κ 0 = ϕ 0 d s δ c s κ 0


(9a) L = 1 - ϕ 0 δ d f c = ϕ 0 δ d s c s .

Equation (9a) shows that L scales both with the geometric mean of df and δ at small melt fractions and with the geometric mean of ds and δ at large melt fractions. Thus, L is a natural length scale associated with thermal equilibrium of fluid-filled pores. The above scaling laws for S justify using the term ϕ0(1−ϕ0) in the scaling length L. It should be noted that we introduce and understand ds as the average distance between melt-filled pores or channels, which can be considerably larger than the grain size. Then both δ and ds and thus L can reach some considerable fraction of the system dimension.

We assume that the fluid and solid phases have the same densities and thermal properties (but relax this assumption later in Sect. 5.1.3):

(13) c p , f = c p , s = c p , 0 , ρ f = ρ s = ρ 0 , κ f = κ s = λ eff c p , 0 ρ 0 = κ 0 .

From Eqs. (4), (6), and (7) we get the non-dimensional energy conservation equations for the fluid and solid phases, respectively:


From these equations we notice that the thermal evolution of the two-phase system is controlled by two non-dimensional numbers: the scaling porosity ϕ0 and the Péclet number Pe defined as

(16) Pe = v f 0 L κ 0 .

This number has already proven to be of high significance for determining whether thermal non-equilibrium is present or not (Schmeling et al., 2018), and the highest Pe corresponds to the largest temperature difference between fluid and matrix. In the following we drop the primes, keeping all equations non-dimensional, if not indicated otherwise.

In the following we consider a homogeneous two-phase matrix–fluid system in 1D with a porosity constant in space and time, i.e., ϕ=ϕ0. We assume a constant fluid velocity which will be expressed in terms of Pe, thus we choose the non-dimensional velocity vf=1. This simplifies Eqs. (14) and (15) to

(17) T f t + Pe T f z = 2 T f z 2 - 1 - ϕ 0 T f - T s


(18) T s t = 2 T s z 2 + ϕ 0 T f - T s ,

respectively. As we are interested in the evolution of the non-equilibrium temperature difference between the solid and fluid, subtraction of Eq. (18) from Eq. (17) gives

(19) T f - T s t - 2 T f - T s z 2 + Pe T f z + T f - T s = 0 ,

which is equivalent to

(20) T f - T s t - 2 T f - T s z 2 + Pe T f - T s z + T f - T s = - Pe T s z .

Note that while the temperatures Tf and Ts explicitly depend on two non-dimensional numbers Pe and ϕ0, the temporal evolution of the temperature difference (TfTs) explicitly depends only on Pe. However, implicitly it is still a function of ϕ0 because Ts on the right-hand side of Eq. (20) depends on ϕ0 via Eq. (18). Only for cases or stages with Ts independent of ϕ0 as proposed in Sect. 4 is the temperature difference (TfTs) a function of only one non-dimensional parameter, Pe, and not of ϕ0.

2.3 Model setup

The fluid and solid heat conservation equations are solved in a 1D domain of height H. Other geometries could also be easily explored but are not considered here, since we focus on studying the relative control of the scaling parameters on thermal non-equilibrium evolution. At time t<0, both solid and fluid are at rest, in equilibrium. Both initial temperatures decrease linearly from 1 to 0 with z; therefore a constant temperature gradient of -G=-1/H is present in both phases (see Fig. 1). As a boundary condition both phase temperatures are set to 1 (non-dimensional temperature) at z=0. At z=H a constant thermal gradient condition T/z=-1/H (non-dimensional) is imposed for both phases. At z=0 the advective flux is fixed by the constant temperature condition (i.e., it is equal to Peϕ0) while at z=H it evolves freely with the fluid temperature (i.e., it is given by TfPeϕ0, all non-dimensional). This top boundary condition needs some justification: the hyperbolic partial differential Eqs. (17) or (18) require two well-defined boundary conditions each, Dirichlet (fixed temperature), Neumann (fixed thermal gradient), Robin (linear combination of Neumann and Dirichlet), or Cauchy (fixed temperature and thermal gradient). Applying the Dirichlet condition at the bottom leaves either a Dirichlet, a Neumann, or a Robin condition to specify for the top. Different combinations of these boundary conditions can be applied separately for the fluid and the solid. For example, for the fluid a Neumann condition with zero or a small temperature gradient may be reasonable, while for the solid one may consider a Robin boundary condition mimicking a thick conductive lid with an internal constant temperature gradient and a fixed surface temperature. However, temporal changes of the temperature at the top (for our system, equal to the bottom of the lid) would lead to unphysical variations in the constant slope of the temperature gradient within the imagined lid. This is because the temperature in the lid can only vary on the diffusive timescale of the lid, which is much longer than all timescales in our model as long as the lid is thicker than H. In fact, in an open outflow situation like our system the evolution of the temperature, the thermal gradient, and the total (advective plus conductive) heat flux are not known a priori, but depend on the evolution within the system. In the early stage of model evolution, both the solid and fluid have a thermal gradient inherited from the initial condition, which is advected upwards in the fluid. Thus it seems most appropriate to use the Neumann condition as a boundary condition for both the solid and fluid. Only at later stages does this boundary condition impose artifacts in the temperature field close to the top boundary. The limitations of this top boundary condition are tested and discussed in Sect. 5.1.2.

Figure 1Initial and boundary conditions.


This model setup adds a third non-dimensional scaling parameter to the system, namely G=1/H. It defines the initial non-dimensional temperature gradient or conductive heat flux, positive for a flux directed upwards. To summarize, the temperatures depend on the non-dimensional parameters Pe, ϕ0, and G.

2.4 Numerical scheme

The equations are solved by a MATLAB (MATLAB R2021b) code using a finite difference scheme central in space for the conduction terms, upwind for the advection term, and explicit in time. The spatial resolution is dz=0.1 or, for a few cases in Fig. 3 below, min0.1,H/100 for H<10. The time step was chosen as dt=14mindz/Pe,dz2, i.e., taking the minimum of the Courant or diffusion criterion. Tests with higher spatial and temporal resolution have been carried out and did not change the results visibly. The global heat balance has also been checked: the maximum relative heat balance error can be defined as δq=qtotz=0-qtotz=H-HTmeantqtotz=0+qtotz=H/2max , where qtot is the total non-dimensional vertical heat flux (conductive and advective) and Tmean is the mean temperature of the model. δq has an error on the order of 1 (due to the upwind scheme) with respect to the grid size dz; i.e., it is approximately equal to const ⋅ dz, where the constant is of the order of 0.2 (i.e., 2 %) for high Péclet numbers and drops to 0.1 (1 %) or 0.01 (0.1 %) for Pe≅1 or smaller, respectively.

3 Numerical model results

First, some example numerical results are shown in Fig. 2 to understand the physics and the typical behavior.

Figure 2Typical model evolution for Pe=1, two different melt fractions ϕ, and two different non-dimensional temperature gradients G (i.e., heights H). (a) Model 1 is with G=0.1 (H=10) and ϕ=0.1. Red and blue curves show the fluid and solid temperatures at different non-dimensional times t as indicated by the legend, respectively. Initial temperatures are almost identical to the t=0.5 curves. Steady state is reached at about t=100; the curves of the last two times plot on each other. (b) Model 2 with G=0.01 (H=100), otherwise as in (a). Steady state is not fully reached. (c) Conductive (blue and red curves) and advective (magenta curve) heat fluxes through the solid and fluid, respectively, of model 1. Line styles indicate the same times as in (a). (d) Same as (c) but for model 2. (e) Temporal evolution of fluid and solid temperatures, Tf (red) and Ts (blue), respectively, at the top of model 2 with ϕ=0.1 and model 3 with ϕ=0.2. G=0.01 (H=100) for both models. (f) Evolution of fluid–solid temperature difference (TfTs) at different distances z in model 2 (ϕ=0.1, solid curves) and in model 3 (ϕ=0.2, dashed curves). (g) Zoomed-in early temporal evolution of solid and fluid temperatures of models 2 and 3 shown in (e). (h) Zoomed-in early temporal evolution of temperature difference of models 2 and 3 shown in (f).


3.1 Evolution of temperatures and thermal non-equilibrium with time

Three different models have been run, all with Pe=1 and the following other parameters:

  • model 1 – G=0.1 (H=10), ϕ=0.1,

  • model 2 – G=0.01 (H=100), ϕ=0.1, and

  • model 3 – G=0.01 (H=100), ϕ=0.2.

Figure 2a and b show Tf and Ts as functions of z at different times as indicated for two initial temperature gradients, G=0.1 (H=10) and G=0.01 (H=100). Figure 2c shows the different contributions to the depth-dependent conductive and advective heat fluxes through the solid and fluid phases, respectively. Figure 2e shows the evolution of Tf and Ts with time at the top of the domain, for the same model 2 as in Fig. 2b and for model 3 with a higher melt fraction ϕ=0.2. Figure 2f shows the evolution of (TfTs) at different distances z of model 2 (ϕ=0.1) and of model 3 (ϕ=0.2). At each depth of the system, the fluid and solid temperatures, as well as the temperature difference and the heat fluxes, evolve following three stages.

Stage 1. During this transient stage the fluid temperature increases faster than the solid temperature (Fig. 2a, b, e, f), and the temperature difference (Fig. 2f, h) increases. During this stage, the fluid temperature increases rapidly at first, and then the temperature increase slows down. The conductive heat fluxes in both solid and fluid decrease rapidly and more slowly later, while the advective heat flux rapidly increases. As for the solid temperature, it first increases slowly and then faster and faster. At t=0, the fluid velocity is suddenly set to non-zero; thus the fluid temperature starts to deviate from equilibrium and increases due to these new conditions. If the solid temperature were maintained constant with time, the fluid temperature would probably reach a steady-state profile, depending on boundary conditions, fluid velocity, and solid temperature. While the fluid temperature increases faster than the solid temperature, the fluid–solid temperature difference, and thus the heat transfer term, increases too, forcing the solid temperature to progressively increase. At the end of stage 1 the maximum temperature difference is approached (Fig. 2h). Because the solid temperature has not risen significantly at that time (at t=4 in the example) compared to the fluid temperature (Fig. 2g), different melt fractions do not affect the temperature differences during this stage (Fig. 2h in which all curves merge into one curve). This observation confirms the expectation from Eq. (20) that the temperature difference does not depend on melt fraction as long as the solid temperature is independent of ϕ, which is the case as long as Ts stays close to its initial profile.

Stage 2. The fluid and the solid temperatures increase at similar rates, constant with time (Fig. 2e), and the temperature difference remains constant and at a maximum at the top (Fig. 2f). Solid–fluid heat transfer is at a maximum during this stage. As Ts is no longer constant in time, different melt fractions lead to different rates of temperature increase (Fig. 2e) and also to different evolutions of (TfTs) (Fig. 2f solid curves compared to dashed curves). At higher melt fraction the heat transfer into the solid increases (see last term in Eq. 18), resulting in a faster increase in the solid temperature whose gradient flattens earlier. Thus, the end of stage 2 is reached earlier (Fig. 2e).

Stage 3. As the fluid temperature rises closer to the Tf value at the bottom, its increase slows down, and heat transfer, and thus temperature difference, decreases. In model 1 (Fig. 2a), steady state is reached while the fluid and solid temperatures are still far from 1. This is due to the influence of boundary conditions, as the heat transferred from the fluid phase to the solid phase is compensated for by the solid phase heat loss at the top of the domain. In model 2 (Fig. 2b), boundary conditions at z=H are applied farther away from the bottom, therefore allowing for a higher increase in temperatures when reaching the steady state.

At each z we observe that the temperature difference first increases rapidly to reach a maximum after a short time (stage 1), hereafter t=4 (Fig. 2h). The resulting amplitude of the temperature difference is identical at the different z positions and for both melt fractions. Then it stays constant at this maximum value (stage 2) and finally decreases (stage 3) (Fig. 2f). The higher in the model, the longer the temperature difference remains at maximum. A higher melt fraction accelerates the decrease in (TfTs).

The absolute maximum temperature difference in space and time does not depend on boundary conditions (see also Sect. 5.1.2 where the influence of boundary conditions is discussed) or on the z position or on the melt fraction and therefore looks to be an interesting observable. It could indeed be useful for getting a first-order estimate of thermal non-equilibrium conditions and possible temperature differences in a magmatic system. In the following sections we study how this maximum temperature difference evolves when varying the parameter Pe.

Comparing the heat fluxes of model 1 (G=0.1) with those of model 2 (G=0.01) shows the importance of heat advection by the fluid phase: in model 1 (Fig. 2c) the conductive contribution through the solid is of the same order of magnitude as the advective contribution by the fluid because the initial temperature gradient G and the porosity φ are the same (=0.1). In model 2 (Fig. 2d), in spite of the same Péclet number, the smaller initial temperature gradient (G=0.01) reduced the conductive contribution with respect to the advective contribution by a factor of about 10, and advection dominates. Furthermore, the two models demonstrate that the conductive heat flux contributions may be important with respect to advective and interphase heat flux for sufficiently large initial thermal gradient (here 0.1), while it has been neglected in several earlier investigations (e.g., Schumann, 1929; Spiga and Spiga, 1981).

3.2 Maximum temperature difference

The maximum temperature difference of a model can be defined as the maximum value reached in space and time (see Fig. 2f). A series of models have been carried out for the two different non-dimensional parameters Pe and G=1/H, and ΔTmax has been determined for each model (Fig. 3). Some first observations can be made.

  • For all Pe values, ΔTmax is proportional to Pe (Fig. 3a) as long as ΔTmax is somewhat smaller than the absolutely possible maximum 1, which is asymptotically approached for high Pe.

  • ΔTmax is proportional to G, i.e., to the non-dimensional temperature gradient for G<0.1.

  • ΔTmax reaches a maximum for large G of the order of 1, i.e., when H reaches 1 or the dimensional H reaches the length scale L.

  • ΔTmax is essentially independent of ϕ as models with different ϕ values almost merge in the same points shown in Fig. 3. This has been verified by running all models of Fig. 3 with melt fractions between 0.1 and 0.9 (not shown).

These observations suggest the existence of several domains in which scaling laws for ΔTmax could be derived, based on the two scaling parameters. In the next section, we propose an analytical derivation of ΔTmax values to obtain scaling laws and confirm the observed proportionalities.

Figure 3Maximum fluid–solid temperature differences TfTs of numerical models (asterisks) with different parameters, plotted (a) as a function of the Péclet number Pe for G=0.1 and ϕ=0.1 and (b) as a function of the initial thermal gradient G for Pe=1 and ϕ=0.1. The solid lines give the analytic solutions. The inset in (a) shows the comparison of the analytic solution Eq. (22) with the different limits derived in Sect. 2.1 and 2.2. The black curve represents the analytic solution, and the colored straight lines show the results in the high or low value limits of Eqs. (26) to (30), respectively. A larger version of the inset is given as Fig. S1 in the Supplement. The analytical solution for the full parameter range PeG is given in Fig. 4.


4 Scaling laws derived from analytical solution

In this section a simplified analytical solution for the z-dependent temperature difference between fluid and solid will be derived. From this solution the maximum temperature differences ΔTmax can be obtained, and scaling laws will be derived.

4.1 Analytical solution of the governing equations

We are interested in an analytical solution of Eq. (20) controlling the non-equilibrium temperature difference (TfTs). We simplify the problem by considering the hypothetical case in which (TfTs) does not change with time and, moreover, in which the thermal gradient in the solid phase is fixed and linear, with Ts/z=-G=-1/H (non-dimensional, with dimensions G=ΔT0/H). Although different from initial and steady-state stages described in the 1D models (Sect. 3.1), this hypothetical case is quite similar to what can be observed at the very beginning of the second stage described in Sect. 3.1 (see Fig. 2f, h). In this second stage, the evolution of Tf and Ts was observed to be quite similar indeed. Besides, at the end of stage 1 (Sect. 3.1), Ts remains close to initial conditions; therefore a fixed linear gradient of slope -G=-1/H is justified. Since the maximum temperature difference between the two phases is observed starting from the end of stage 1 and during stage 2 (Sect. 3.2), it does not seem unreasonable to consider this hypothetical case for deriving the maximum temperature difference. Using these assumptions, Eq. (20) becomes

(21) 2 T f - T s z 2 - Pe T f - T s z - T f - T s = - Pe G .

While in the general case of Eq. (20) the temperature difference implicitly depends on ϕ0, i.e., on the three non-dimensional parameters Pe, ϕ0, and G, Eq. (21) no longer depends on ϕ0 because we replaced Tsϕ0/z with G, which is independent of ϕ0. Equation (21) is a second-order ordinary differential equation for (TfTs) whose solution can be analytically derived as (see Supplement for details)

(22) T f - T s = α e r 1 z + β e r 2 z + Pe G ,

where r1 and r2 are the roots of the associated equation of Eq. (21):

(23) r 1 = 1 2 Pe - Pe 2 + 4 , r 2 = 1 2 Pe + Pe 2 + 4 .

The parameters α and β are constrained by the boundary conditions Tf-Ts=0 at z= 0 and Tf-Tsz=0 at z=H:

(24) α = Pe G r 2 r 1 e r 1 - r 2 / G - r 2 , β = Pe G r 1 r 2 e r 2 - r 1 / G - r 1 .

The third term in Eq. (22) is a particular solution for Eq. (21).

4.2 Comparison with numerical models

From Eq. (22) the maximum value of the depth-dependent temperature difference (TfTs) can be determined. It is found that the maximum is always at z=H. This value will be denoted as ΔTmax and has been calculated for all parameter combinations used for the numerical models. In Fig. 3 these analytical solutions are plotted as solid lines together with the numerical solutions (asterisks). The agreement is very good. For most cases the differences between the numerical and analytical solutions are well below 1 %; only when ΔTmax reaches values of about 0.6 and higher do the differences become >1 %, up to 6 %. This general good agreement is another justification for using the time-independent Eq. (21) to obtain an analytical solution of an intrinsically time-dependent process as long as we are interested only in the maximum value of (TfTs). Other reasons for the observed differences between the analytical and numerical solutions include numerical errors when determining the particular times when maximum temperature differences are reached, especially for the models which are in the regime close to ΔTmax=1 where the ΔTmax(Pe) curves become non-linear (Fig. 3a).

4.3 Scaling laws for temperature differences at certain parameter limits

The analytical solution for ΔTmax fits very well with our model results and therefore looks to be ideal for getting a better understanding of the relative influences of the two controlling parameters Pe and G, described in Sect. 2.2 and 2.3. The Péclet number is already known to be of great importance for thermal equilibrium and non-equilibrium conditions. Inspecting the last term in Eq. (22) we notice that a high Pe and a high initial thermal gradient should favor higher temperature differences. This has been demonstrated in Fig. 3.

Equation (22) is, however, complicated, and the assessment of the relative importance of Pe and G for different possible regimes is limited. In this section, we study the evolution of (TfTs), i.e., also ΔTmax, in a few limiting cases. This enables us to better understand the influence of each parameter and to derive some scaling laws for different regimes.

4.3.1 Limit Pe→0

When Pe tends to 0, we have the condition

(25) Pe 2 .

With this condition Eq. (22) tends to the following limit (see Supplement):

(26) T f - T s = Pe G 1 - M ,


(27) M = cosh ( z ) + cosh 2 G - z 1 + cosh ( 2 / G ) ,

which simplifies for z=H=1/G to

(28) M = 1 cosh ( 1 / G ) .

This is the limit for Pe→0. This limit gives predictions for ΔTmax in very good agreement with Eq. (22) for Pe<1 (having G=0.1) (see inset in Fig. 3a or Fig. S1 in the Supplement). In the limit G→0 and finite Pe<1/G we get the limit for M


Thus, for both small Pe and small G the temperature difference (Eq. 26) can be written as

(29) T f - T s = Pe G 1 - e - z .

Equation (29) confirms the proportionalities observed in Fig. 3, namely ΔTmaxPe (Fig. 3a) and ΔTmaxG (Fig. 1b).

4.3.2 Limit Pe→∞

To obtain the limit of Eq. (22) for Pe→∞, Eq. (22) can be linearized with respect to 4/Pe21. Applying the rule of L'Hôpital, Eq. (22) tends to the following limit:

(30) T f - T s = G z .

For details, see the Supplement. This limit is also the solution of Eq. (21) when neglecting the diffusive and heat transfer terms. As demonstrated in the Supplement, this limit predicts ΔTmax values in very good agreement with Eq. (22) for Pe>100 (Fig. 3a, inset).

4.3.3 Exploring the domains for the maximum temperature difference including all limits

Before exploring the full parameter space, we first give a short overview of expected parameter ranges in magmatic systems.

In natural magmatic systems such as mid-ocean ridges, Pe is expected to evolve from very low values of the order of 10−5 to 10−3 in partially molten regions with distributed porous flow to higher values of the order of 1 or larger at depths where channels have merged and further to very high values of the order of 105 in dike systems (Schmeling et al., 2018).

While the melt fraction does not influence ΔTmax (cf. Eqs. 22, 30), it does influence the long-term temporal behavior because Ts is ϕ0-dependent (see Eq. 20). Therefore, here follow some words about possible melt fractions. As melt flow may occur at very small melt fractions (McKenzie, 2000; Landwehr et al., 2001), large ϕ values are not expected in natural mantle magmatic systems or in dike systems in the crust. Values of channel volume fraction generally remain below a few percent up to tens of percent (in dunite channels up to 10 %–20 %, Kelemen et al., 1997).

To get an idea about the expected order of magnitude of the macroscopic thermal gradient G=1/H of the system, we have to evaluate the scaling length L used to scale the dimensional H. L scales with the geometric mean of the channel width df and the interfacial boundary layer thickness δ (Eqs. 9 with 11). L would evolve non-linearly with the width of melt pathways, which may increase by several orders of magnitude as 3D grain junctions eventually merge to 1D dikes. As will be shown in Sect. 5.3 in more detail the resulting non-dimensional G ranges between an order of 1 and an order of 10−5.

In Fig. 4 we explore ΔTmax variations using the analytical solution Eq. (22), in which ΔTmax depends on Pe and G. Three main regimes can be distinguished.

  • Regime 1. For high Pe values, (TfTs) tends to the relationship described in Eq. (30). The temperature difference increases linearly with distance from the bottom (z=0), reaching ΔTmax=1 at z=H. In the whole region the fluid temperature remains constant and at a maximum of 1 while the solid temperature increases linearly with z from 0 to 1. The proportionality of ΔTmax to G disappears because the maximum value of z is equal to H=1/G.

  • Regime 2. For Pe≪1, or more precisely, for Pe1G represented by the oblique dashed line in Fig. 4, (TfTs) varies with distance from the bottom according to 1-e-z and is proportional to Pe and G. This means that large temperature gradients favor large temperature differences. In this domain, (TfTs) tends to the relationship presented in Eq. (29).

  • Regime 3. For a large initial temperature gradient G close to 1 (small H) and Pe≪1, (TfTs) tends to the relationship proposed in Eq. (26). In this domain, (TfTs) is proportional to Pe but no more to G because M is a function of G gradually canceling the proportionality to G, which is visible in regime 2. The depth dependence is given by (1−M(z)), which at G=1 increases non-linearly from about 0 to 0.4 with increasing z.

Figure 4Main regimes of the maximum fluid–solid temperature differences ΔTmax due to thermal non-equilibrium obtained by the analytical solution (Eq. 22) in the parameter space of the Péclet number Pe and temperature gradient G. The asymptotic limits are indicated by the formulas, and M(z) is given by Eq. (27) with 1-MG,z increasing non-linearly from about 0 to about 0.4 with increasing z for G in the range 0.3 to 3. Regime boundaries are shown as dashed lines. Typical parameter combinations for magmatic settings such as interstitial melts or dikes are indicated by the orange rectangles which extend further to the left, well below log 10G of −3. Note that two slices through this field at G=0.1 and at Pe=1 have already been shown in Fig. 3.


5 Discussion

5.1 Limitations

5.1.1 Comments on the analytic solution

Although the assumptions used to get the analytic solution (Eq. 22) are very specific, they are reasonable considering the conditions in the models when ΔTmax is reached, and it fits the numerical results very well. This is shown in Fig. 5 where for various combinations of Pe and G the time-dependent temperature differences (TfTs) are shown as functions of depth together with the analytical solutions using Eq. (22). For all examples the position of the maximum temperature differences lies at z=H. A major simplification used in Eq. (21) was time independence. Obviously, the resulting analytical solutions represent stage 2, which is quasi steady state in contrast to stage 1 when the temperature difference builds up and stage 3 when the long-term behavior is approached. We emphasize that this analytical solution is a very good approximation of the depth-dependent temporal maximum temperature difference that can be reached in such porous systems.

Figure 5Comparison of depth- and time-dependent numerical solutions with the time-independent analytical solutions for different parameters Pe and G as indicated in the panel titles. In each panel the curves show (TfTs) profiles for progressive times, and the colors are cyclically varied with time from blue to yellow, starting with blue (bold curve). The bold red dashed curve shows the analytical solution Eq. (22), which represents a very good estimate of the depth-dependent temporal maximum of the temperature difference. In each panel the first five curves are plotted at time increments of 0.5 (0.025 for Pe=100) and the later curves with 5 (1 for Pe=100). The total non-dimensional times of each panel are 100 (500 for G=0.01). Steady state is reached in the models with G=0.1. The model with G=0.01 would need t= 10 000 to reach steady state. The melt fraction was chosen as ϕ0=0.1.


5.1.2 Boundary conditions at top and initial conditions

The boundary conditions we chose at the top (z=H) are suitable for cases with little temperature evolution (regimes 2 and 3, low Pe) and for early stages for regime 1 but might be inappropriate for high temperature increases (high Pe – regime 1) at later stages (see Sect. 4.3.4). In order to quantify the influence of this choice of boundary conditions on our results, we compared the evolution of (TfTs) profiles for three Péclet numbers and two values of G, using four different boundary conditions at the top (Fig. 6).

  • Constant thermal gradient equal to the initial thermal gradient in the solid and fluid phases (Neumann condition). This was the boundary condition used in the models.

  • The thermal gradient is set to 0 at the top (Neumann condition).

  • Both fluid and solid temperatures are set to 0 at the top (Dirichlet condition).

  • Temperature at the top is numerically calculated from the full Eqs. (17) and (18) using one-sided (upwind) positions for the first and second derivatives (open boundary).

Mathematically, the open boundary condition is not a rigorous boundary condition because both the temperature and temperature gradient intrinsically depend on the temperature evolution within the model. Therefore, it cannot be applied to the analytical solution of Sect. 4.1. Numerically it works well for our system without producing instabilities or oscillations. Comparing the top and bottom rows of Fig. 6, the constant temperature gradient condition produces quite similar results as the open boundary condition for all Pe and G values tested during the first and second stages of temporal evolution (see Sect. 3.1). The agreement becomes worse for stage 3 when approaching steady state at large Pe values. Comparing the other two boundary conditions (second and third rows of Fig. 6) with the constant gradient condition (top row) shows that the effect of the top boundary during stages 1 and 2 is still small sufficiently far away from the top. Only for the small Pe case (left column of Fig. 6) do the zero gradient and zero temperature conditions strongly affect the upper half of the domain by diffusion. Yet the maximum temperature difference of the constant gradient case is nearly reached by the other two boundary conditions further within the domain, not at the top. The special case of high Pe and small G with zero temperature boundary condition (third row, fourth column in Fig. 6) shows a strong build-up of TfTs close to the top when approaching the steady state. This stems from the large local temperature gradient built up near the top as a result of transforming the difference in advective heat in- and output PeTinflux-PeToutflux=Pe into a high conductive outflux T/z at the top. It is unlikely that such situations occur in natural systems.

Figure 6Temporal evolution of vertical profiles of (TfTs) for models with different Péclet numbers and different initial temperature gradients G. In each panel the curves show (TfTs) profiles for progressive times, and the colors are cyclically varied with time from blue to yellow, starting with blue (bold curve). The first five curves of the Pe<100 (respectively Pe=100) models were taken with time increments of 1 (respectively 0.1) and the later curves with 10 (respectively 1). The total time was 100 in all models with G=0.1 and 500 in the models with G=0.01. Steady state is reached in the models with G=0.1. The models with G=0.01 would need t= 10 000 to reach steady state. In each row the top boundary conditions are assumed as indicated at the left.


We have also tested an open boundary condition for the fluid and a Robin boundary condition for the solid imagining a lid on top of our model with a constant temperature gradient and a fixed surface temperature. Choosing the surface temperature in such a way that the initial thermal gradient within our system and within the lid are identical, this boundary condition adds a new non-dimensional parameter, the lid thickness Hlid. Tests show that lid thicknesses larger than H give results in general agreement with the constant gradient or open boundary models of Fig. 6. Localized differences with up to 30 % higher (TfTs) values occur near the top, but they disappear for Hlid≥10 H. Lid thicknesses smaller than H force the solid temperature at the top to values closer to 0, while the fluid temperature remains high. This results in significantly higher (TfTs) values than in Fig. 6 (top row) close to the top. However, typical natural magmatic systems are on the Hlid>H side, suggesting that our constant gradient boundary condition is a good approximation.

In summary, the influence of boundary conditions on fluid and solid temperature evolution depends mostly on the domain size H and on the value of Pe. The larger these two parameters, the less important the influence of boundary conditions within almost the whole model domain. If one is interested in the maximum value of TfTs in space and time, the tests show that this value can safely be picked at z=H when using the constant temperature gradient boundary condition.

As an initial condition we used a linear temperature profile and initial equilibrium between solid and fluid. A non-linear initial temperature profile between Tf=Ts=1 at the bottom and Tf=Ts=0 at the top would have spatially varying temperature gradients with sections with gradients larger than those assumed in our model. As the temperature gradient strongly influences thermal non-equilibrium (see e.g., Eq. 22, which explicitly contains the temperature gradient G), the above results are expected to be different, and a stronger thermal non-equilibrium is expected in regions with higher gradients. Schmeling et al. (2018) used a step function with Tf=Ts=1 at z=0 and Tf=Ts=0 at z>0 as the initial condition, i.e., an extremely non-linear profile near z=0. Assuming this initial temperature profile, Fig. 7 shows the temporal behavior of the temperature difference for selected parameter combinations, equal to the parameters used in Fig. 5. The analytical solutions for the time-independent cases (Eq. 22) are also shown. As expected, at early stages the temperature differences are significantly larger than given by the analytical solutions by a factor of 2 or more shortly after the onset of the evolution. At later stages (stage 2 or 3) the time-dependent solutions approach or pass through the analytical solutions. Thus, we may state that the analytical solutions depicted in the regime diagram in Fig. 4 represent lower bounds of thermal non-equilibrium compared to settings with non-linear initial temperature profiles.

Figure 7Time- and depth-dependent numerical solutions (thin curves) as in Fig. 5 but for step-function initial conditions: Tf=Ts=1 at z=0 and Tf=Ts=0 at z>0 at t=0. The bold dashed red curves are the time-independent analytical solutions as in Fig. 5. In each panel the curves show (TfTs) profiles for progressive times, and the colors are cyclically varied with time from blue to yellow, starting with blue (bold curve). In each panel the first five curves (and later curves, respectively) are plotted at time increments of (a) 0.5 (5), (b) 1 (10), (c) 0.5 (5), and (d) 0.025 (1). The total non-dimensional times of each panel are 100 (500 for G=0.01). Steady state is reached in the models with G=0.1. The model with G=0.01 would need t= 10 000 to reach steady state. As porosity ϕ=0.1 is assumed.


5.1.3 Different densities and thermal properties of the two phases

While for simplicity we used equal physical properties for the fluid and solid, in many circumstances they might be significantly different. Equal properties are good approximations for magmatic systems where differences of density and thermal parameters are small (order of 10 %), whereas porous flows of water or gases through rocks or other technical settings may be characterized by larger differences. Allowing for different material properties adds four new parameters, namely the ratio of diffusivities, the ratio of densities, the ratio of heat capacities, and a new effective thermal conductivity λeff for the interface between the two phases with different properties. To evaluate how many new non-dimensional numbers are introduced, we non-dimensionalize the equations assuming different material properties for the two phases. We use the fluid properties as scaling quantities and assume that they are independent of temperature, pressure, and depth. We modify the scaling length and Péclet number by defining L̃=Lλeff with λeff=λeff/λf and Pẽ=Peλeff. With this scaling, Eqs. (14) and (15) turn into (for clarity, primes indicate non-dimensional quantities)

(31) ϕ T f t + Pe ̃ v T f = ( ϕ T f ) - ϕ 0 1 - ϕ 0 ( T f - T s )


(32) ( 1 - ϕ ) T s t = κ s ρ s c p , s ( ( 1 - ϕ ) T s ) + ϕ 0 1 - ϕ 0 1 ρ s c p , s ( T f - T s ) .

Inspection of these equations shows that two more non-dimensional numbers are introduced: the ratio of diffusivities κs and the ratio of the products' density and heat capacity, ρscp,s.

As Eqs. (31) and (32) cannot be merged into one time-independent ordinary differential equation for (TfTs) as in Sect. 4.1, we numerically tested some cases with Pẽ=Pe=1 and λeff=1 in which the diffusivity ratio and the ratio of ρscp,s were varied between 0.1 and 10 (see Fig. 8). The results show that for the fixed combination of Pe=1 and λeff=1 the magnitude of thermal non-equilibrium remains on the same order of magnitude O(0.1) as for equal properties (Fig. 8). However, the time dependence is significantly affected: for a high ratio of κs=10 (i.e., the solid is strongly conducting) the solid temperature profile remains close to the constant initial gradient, and the temperature difference rapidly converges to a steady state similar to the analytical solution depicted in Fig. 5a. In contrast, for a low κs=0.1 the solid temperature departs more strongly from the initial linear gradient, and the solid–fluid temperature difference slowly drops with time in the long term. Varying the potential to store heat in the solid, i.e., ρscp,s, Fig. 8e and f show that a high value slows down the long-term time-dependent variations, while a small value leads to rapid long-term temporal variations in (TfTs) and faster convergence to the steady state, which is similar to the case with equal properties.

Figure 8Time- and depth-dependent profiles of the fluid–solid temperature differences as in Fig. 5. (a) Reference models (as in Fig. 5a) with Pe=1, G=0.1, ϕ=0.1, and equal fluid-to-solid properties. (b–f) Profiles as in (a) but with solid-to-fluid property ratios as indicated in the titles of each panel, and λeff=1. The properties in (b) are typical for water in sedimentary rocks. In every panel but (b) the first five curves were taken with time increments of 0.5 and the later curves with 5. In (b) the first 5 curves were taken with time increments of 0.4875 and the later curves with 4.875. The total time was 100 in all models. Steady state is reached in all models.


It is interesting to apply the results for different physical properties to a geologically relevant setting, namely water flowing through sedimentary rocks. Given that the high heat capacity of water is about 3 times larger than that of rock, and the density is almost 3 times less, the product ρscp,s is about 0.78, i.e., of an order of 1. However, the thermal diffusivity of water is significantly smaller than that of rock, typically by a factor of 16; i.e., κs is about 16. We tested a few cases (Fig. 9) with Péclet numbers and initial thermal gradients G (i.e., inverse model heights) (assuming for simplicity λeff=1) equal to the cases depicted in Fig. 5. The time-dependent profiles behave similarly to those in Fig. 5, with very similar maxima of the temperature differences (red dashed curves in Fig. 5) relevant for stage 2. The only important difference is that the water–sedimentary rock case more rapidly approaches the late steady states of stage 3, and these stages are closer to the maximum red dashed curves. These results suggest that the absolute values of maximum thermal non-equilibrium temperature differences shown in the regime diagram Fig. 4 are also applicable to a water–sedimentary rock system.

Figure 9Time- and depth-dependent profiles of the fluid–solid temperature differences as in Fig. 5, but for fluid-to-solid property ratios typical for water flowing through sedimentary rocks, i.e., ρscp,s=0.78, κs=16, and λeff=1. Pe and G have been chosen as indicated in the panel titles (as in Fig. 5), and ϕ=0.1 was assumed. In each panel the curves show (TfTs) profiles for progressive times, and the colors are cyclically varied with time from blue to yellow, starting with blue (bold curve). The first five curves were taken with time increments of 0.4875 and the later curves with 4.875. The total time was 100 in all models with G=0.1 and 200 in the models with G=0.01. Steady state is reached in the models with G=0.1. The model with G=0.01 would need t= 10 000 to reach steady state.


5.2 Timescales

It is interesting to evaluate the timescales for reaching the maximum non-equilibrium temperature differences and the steady state. For every numerical model, we recorded the time needed to reach 90 % of the maximum temperature differences between fluid and solid, t90 %, and the time needed to reach steady state, tsteady. The latter has been determined as the time at which the maximum difference between (Tf(z)−Ts(z)) curves at two subsequent time steps becomes less than 10-8ΔTmax. These times can be compared with different timescales that may characterize the evolution of temperatures in the models. These timescales can be based on advection over a characteristic distance dchar, giving tadvd=dchar/vf0, or on diffusion over the characteristic distance, giving tdiffd=dchar2/κ0. We tested these timescales with the two natural length scales of the models. The first is the scaling length L (equal to 1 non-dimensional), representing essentially the geometric mean of the channel width of the pores, df, and the interfacial boundary layer thickness δ. The second is the model height H. Grouping the models depending on the regime they belong to (see Sect. 4.3.4 and Fig. 4), we plotted the recorded times t90 % and tsteady versus the characteristic timescales mentioned above. Good agreement with the characteristic timescales is indicated by observed times fitting to the dashed x=y lines (Fig. 10).

Figure 10For evaluating timescales the numerically determined times of models with various parameters Pe and G representing the three different regimes 1, 2, and 3 (different symbols) are plotted against characteristic scaling times. (a) Times for reaching 90 % of the maximum temperature difference ΔTmax are plotted against either the advective timescale tadvH based on model height H for regime 1 models or against the scaling time t0 for regime 2 models, or against the diffusive timescale tdiffH based on the model height H. (b) Times for reaching steady states are plotted against the characteristic diffusive timescales, tdiffH, based on model height H for all three regimes. Models close to the dashed line (y=x) are in best agreement with the characteristic times. In (a) regime 2 times are taken dimensionally by multiplying the observed times and the non-dimensional scaling time t0=1 by some arbitrary dimensional times t0.


In regime 1 (high Pe), t90 % is proportional to tadvH (Fig. 10a, blue circles). In this regime the high value of Pe makes the fluid temperature increase fast. It reaches its maximum value during the time under which significant fluid–solid heat transfer builds up and the solid temperature is still low. This corresponds to the time for traveling the full distance H. During stages 2 and 3 the solid temperature increases and the temperature difference decreases before steady state is reached. The time for reaching steady state (Fig. 10b, circles) varies roughly linearly with tsteadytdiffH. For most cases it is controlled by diffusion through the solid over distances of the order of H. The case with large H (circle in Fig. 10b below dashed line) apparently reaches the steady state earlier, but still later than on a corresponding advective timescale based on H (not shown). Inspecting this model shows that during stages 2 and 3 the high Pe number facilitates approaching thermal equilibrium rapidly within large parts of the model and reducing the effective length scale (and characteristic timescale) over which non-equilibrium is still present.

In regime 2 (low Pe and G<0.1, i.e., H>10) the time for reaching ΔTmax is controlled by interfacial heat transfer (Fig. 10a, red asterisks) on the length scale L, resulting in t90 % proportional to t0. The time for reaching steady state is controlled by the diffusion timescale across the height of the system (Fig. 10b).

In regime 3 (low Pe and high G (small H)), time for reaching ΔTmax is similar to or shorter than the diffusion time based on the model height H (Fig. 10a, black crosses). The flattening of the curve indicates that non-equilibrium is reached faster for some models because Pe reaches the order of 1 and the advective timescale starts to take over. The time for reaching steady state (Fig. 10b, crosses) varies linearly with tsteadytdiffH. Clearly, it is also controlled by diffusion.

5.3 Applications to magmatic systems

We now test the possible occurrence of thermal non-equilibrium in natural magmatic systems based on the suggested controlling non-dimensional parameters, namely the Péclet number Pe, the initial thermal gradient G (=1/H), and the melt fraction ϕ. Typical stages of melt flow for mid-ocean ridges include three stages:

  • a.

    partially molten regions with interstitial melts sitting at grain corners, grain edges, or grain faces with low (0.0001 %–6 %) melt fractions (see, e.g., the discussion in Schmeling, 2006);

  • b.

    merging melt channel or vein systems with high-porosity (> 10 %–20 %) channels identified as dunite channels after complete melt extraction (Kelemen et al., 1997);

  • c.

    propagating dikes or other volcanic conduits.

Let us assume typical overall melt fractions of 1 % to 20 % for stages (b) and (c). Schmeling et al. (2018) discussed possible Péclet numbers for such systems based on a Darcy-flow-related Péclet number:

(33) Pe D = v D d s κ 0 .

As we preferably use the melt pore dimension df in our scalings (Eqs. 9a and 10a), we need to relate it to the solid phase dimension ds by using

(34) d s = d f g ϕ , g = 1 - ϕ melt channels ϕ 1 - ϕ melt tubes .

Using Eqs. (34), (9a), and (16), we arrive at the Péclet number used here.

(35) Pe = Pe D 1 g c 1 - ϕ δ d f

Schmeling et al. (2018) reviewed and estimated typical pore or channel spacings ds of 10−3–10−2 m for stage (a), 0.1 m for early stage (b) increasing to 1–100 m for late stage (b), and 100–300 m for stage (c) (dikes). Arguing for typical geometries, spreading rates, and melt extraction rates, Schmeling et al. (2018) estimated the Darcy velocity to range between 10−10 and 10−9 m s−1. With these parameters PeD numbers for the three stages can be estimated for the three stages as

  • a.

    10−7 to 10−5,

  • b.

    10−5 to 10−4 at depths where channel distances are of the order of 0.1 m and 10−4 to 0.1 at shallower depths where the channel distances have increased to the order of 1 to 100 m,

  • c.

    >105 for the dike stage.

To estimate Péclet numbers as defined here (Eq. 35), typical interfacial thermal boundary layer thicknesses δ are needed. As the thermal interfacial heat exchange intrinsically is time-dependent, a good estimate is δ=cthκ0t (in dimensional form), where cth is a constant for a thermal boundary layer, equal to 2.32 for a cooling half space (Turcotte and Schubert, 2014). Assuming that the characteristic time can be expressed by the (dimensional) fluid velocity v0 and system height H, i.e., by t=H/v0=Hϕ/vD, we may express vD in terms of the Péclet number PeD. With the resulting t and δ we arrive at the following Péclet number (H and df are dimensional or non-dimensional):

(36) Pe = Pe D 3 / 4 c th c g - 3 / 4 H d f 1 / 4 1 - ϕ .

For mid-ocean ridge settings we assume H to represent the transition region between the lithosphere and asthenosphere with a thickness of the order of 1 to 10 km, and we use Eq. (34) to insert typical df values. They increase from 10−4 m for interstitial melts (stage a) to 10−3 to 10−1 m for the channeling stage (b) (see Schmeling et al., 2018) to >10 m for the dike stage (c). The resulting Péclet number (Eq. 36) is of the order of 10−3 to 0.5 for stage a, 10−2 during the early stage (b) and 10−2 to 1 during the later stage (b) appropriate for dunite systems, and 104 to 107 for the dike stage (c). To estimate typical non-dimensional thermal gradients G (or layer thickness H), the above estimate for δ and df can be inserted into the scaling length L (Eq. 9a) to arrive at a non-dimensional G=1/H.

(37) G = H d s - 3 / 4 g - 1 / 2 ϕ 3 / 4 c th c Pe D - 1 / 4 1 - ϕ

With the derived estimates for the three stages, G is of the order of 10−6 to 10−2.7 for stage (a), 10−4–10−2.5 increasing to 10−4–0.6 for stage (b), and 10−5–10−2 for the dike stage (c). These resulting stages for Pe and G are indicated in the regime diagram (Fig. 4). All three stages extend far into the domain, with G values smaller than 0.001. Thermal non-equilibrium of the three stages can be summarized as follows.

  • a.

    Interstitial melts are at full thermal equilibrium.

  • b.

    Channeling and veining may result in moderate thermal non-equilibrium at sufficiently high thermal gradients.

  • c.

    After transition to diking full thermal non-equilibrium is predicted.

A similar exercise can be done for continental magmatic systems. We skip such an explicit evaluation here but note that silicic melt viscosities are typically higher than those of basaltic melts at mid-ocean ridges. Thus, Péclet numbers are expected to be smaller, but non-dimensional thermal gradients (Eq. 37) might be larger, resulting in a downward and rightward shift of the natural stages indicated in Fig. 4.

To make our scaling laws and timescales for reaching maximum thermal non-equilibrium more accessible, it is worth writing them in dimensional form. First, to estimate the Péclet number of a natural system, combining Eqs. (9) and (16) gives

(38) Pe = v f 0 κ 0 ϕ 0 1 - ϕ 0 δ S ,

indicating that for very small or very large melt fractions Pe becomes very small. One may use Eqs. (11) or (12) to write Pe also in terms of pore or solid (grain or channel spacing) dimensions df or ds, respectively. The scaling laws and characteristic timescales for the three regimes we found (Fig. 4) are in dimensional form.

  • Regime 1. For large Pe values the maximum non-equilibrium temperature difference is simply equal to the imposed temperature difference, ΔTmax=ΔT0, and the characteristic time to reach maximum non-equilibrium is simply tchar=H/vf0, i.e., the total time of a fluid particle for passing through the system.

  • Regimes 2 and 3. For small Péclet numbers (Pe<HSϕ01-ϕ0δ) the maximum temperature difference scales like

    (39) Δ T max = G v f 0 ϕ 0 1 - ϕ 0 δ κ 0 S = v f 0 ϕ 0 1 - ϕ 0 δ H κ 0 S Δ T 0 ,

    and the characteristic time for reaching this non-equilibrium scales with t0, i.e.,

    (40) t char = ϕ 0 1 - ϕ 0 δ κ 0 S .

    These relations can easily be used to assess the potential of thermal non-equilibrium in systems of fluid flow through solids with given geometrical properties and fluid fractions.

In the above discussion we used the terms moderate thermal non-equilibrium, which we may identify with ΔTmax of a few percent to, say, 30 %, while full thermal non-equilibrium includes higher values up to 100 %. To translate this into dimensional ΔTmax values, what are typical ΔT0 values for mid-ocean ridges? In our example we defined Has the thickness of the transition region between lithosphere and asthenosphere. Such a transition zone may be defined by the depth region bounded by the asthenospheric temperature Tasth and, say, 0.8Tasth, i.e., ΔT0=0.2Tasth200 K. Thus moderate non-equilibrium may be represented by excess temperatures of the melt with respect to the solid between, say, 6 and 60 K, while full thermal non-equilibrium suggests the full ΔT0 range of 60 to 200 K or even higher in the case of dikes extending up into the lithosphere with temperatures below 0.8 Tasth. These typical temperature estimates may have some implications for whether the solid rock will melt or the melt will freeze. However, this discussion is beyond the scope of this paper.

6 Conclusions

In conclusion, we showed that in magmatic systems characterized by two-phase flows of melts with respect to solid, thermal non-equilibrium between melt and solid may arise and becomes important under certain conditions. The main conclusions are summarized as follows.

From non-dimensionalization of the governing equations, three non-dimensional numbers can be identified controlling thermal non-equilibrium: the Péclet number Pe, the melt porosity ϕ, and the initial non-dimensional temperature gradient G in the system. The maximum possible non-equilibrium solid–fluid temperature difference ΔTmax is controlled only by two non-dimensional numbers: Pe and G. Both numerical and analytical solutions show that in a PeG parameter space three regimes can be identified.

  • In regime 1 (high Pe(>1/G)) strong thermal non-equilibrium develops independently of Pe, and a non-dimensional scaling law Tf-Ts=Gz has been derived.

  • In regime 2 (low Pe(<1/G) and low G(<0.3)) non-equilibrium decreases proportionally to decreasing Pe and G, and the non-dimensional scaling law reads Tf-Ts=PeG1-e-z.

  • In regime 3 (low Pe(<1) and G of the order of 1) non-equilibrium scales with Pe and G and is depth-dependent. The scaling law is Tf-Ts=PeG1-Mz, where M(z) depends on G.

Further conclusions include the following.

  • The timescales for reaching thermal non-equilibrium scale with the advective timescale in the high-Pe regime and with the interfacial diffusion time in the other two low-Pe-number regimes.

  • Applying the results to natural magmatic systems such as mid-ocean ridges can be done by estimating appropriate orders of Pe and G. Plotting such typical ranges in the PeG regime diagram reveals that (a) interstitial melt flow is in thermal equilibrium, (b) melt channeling as revealed by dunite channels, for example, may reach moderate thermal non-equilibrium, and (c) the dike regime is at full thermal non-equilibrium.

  • In the studied setup G was constant, leading to conservative estimates of thermal non-equilibrium. Any other depth-dependent initial temperature distributions generate higher non-equilibrium than reported here.

  • The derived scaling laws for thermal non-equilibrium are valid for equal solid and fluid properties. Assuming different properties such as for a water–sandstone system results in similar maximum non-equilibrium temperature differences, but in significantly different time evolutions.

While for simplicity the presented approach has been done essentially for constant model parameters, it can easily be extended to vertically varying parameters. Thus, tools are provided for evaluating the transition from thermal equilibrium to non-equilibrium for anastomosing systems (Hart, 1993).

Code availability

The MATLAB code is listed in the Supplement and is available on request.

Data availability

There are no underlying research data other than the graphs given in the figures. These graphs can be retrieved by the MATLAB programs given in the Supplement.


The supplement related to this article is available online at:

Author contributions

LC formulated the scientific question and the theory of the first version and wrote a first draft including model results. HS did the final writing up, improved the scaling, and reran all models.

Competing interests

The contact author has declared that neither they nor their co-author has any competing interests.


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


We gratefully acknowledge the excellent reviews by John Rudge and Cian Wilson, who stimulated us into significantly improving the scaling.

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. 403710316).

Review statement

This paper was edited by Juliane Dannberg and reviewed by John Rudge, Cian Wilson, and Samuel Butler.


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

Amiri, A. and Vafai, K.: Analysis of Dispersion Effects and Non-Thermal Equilibrium, Non-Darcian, Variable Porosity In-compressible Flow Through Porous Media, Int. J. Heat Mass Tran., 37, 939–954, 1994. 

Becker, K. and Davis, E.: On situ determinations of the permeability of the igneaous oceanic crust, in: Hydrogeology of the Oceanic Lithosphere, 311–336, edited by: Davis, E. and Elderfield, H., CambridgeUniv. Press, ISBN 0 521 81929, 2004. 

Bruce, P. M. and Huppert, H. E.: Solidification and melting along dykes by the laminar flow of basaltic magma, in: Magma transport and storage, edited by: Ryan, M. P., Wiley, Chichester, 87–101, ISBN 0 471 92766 X, 1990. 

Davis, E. E., Chapman, D. S., Wang, K., Villinger, H., Fisher, A. T., Robinson, S. W., Grigel, J., Pribnow, D., Stein, J., and Becker, K.: Regional heat flow variations across the sedimented Juan de Fuca ridge eastern flank: constraints on lithospheric cooling and lateral hydrothermal heat transport, J. Geophys. Res., 104, 17675–17688, 1999. 

de Lemos, M. J. S.: Thermal non-equilibrium in heterogeneous media, Springer Science+Business Media, Inc.,, 2016. 

Furbish, D. J.: Fluid Physics in geology, Oxford University Press, New York, 476 pp., ISBN 0-19507701-6, 1997. 

Harris, R. N. and Chapman, D. S.: Deep seated oceanic heat flow, heat deficits and hydrothermal circulation, in: Hydrogeology of the Oceanic Lithosphere, 311–336, edited by: Davis, E. and Elderfield, H., CambridgeUniv. Press, ISBN 0 521 81929, 2004. 

Hart, S. R.: Equilibration during mantle melting: a fractal tree model, P. Natl. Acad. Sci. USA, 90, 11914–11918, 1993. 

Kelemen, P. B., Whitehead, J. A., Aharonov, E., and Jordahl, K. A.: Experiments on flow focusing in soluble porous media, with applications to melt extraction from the mantle, J. Geophys. Res., 100, 475–496,, 1995. 

Kelemen, P. B., Hirth, G., Shimizu, N., Spiegelman, M., and Dick, H. J. B.: A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges, Philos. T. R. Soc. S.-A, 355, 283–318, 1997. 

Keller, T., May, D. A., and Kaus, B. J. P.: Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust, Geophys. J. Int., 195, 1406–1442, 2013. 

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–3933,, 1994. 

Landwehr, D., Blundy, J., Chamorro-Perez, E. M., Hill, E., and Wood, B.: U-series disequilibria generated by partial melting of spinel lherzolite, Earth Planet. Sc. Lett., 188, 329–348, 2001. 

Lister, J. R. and Kerr, R. C.: Fluid-mechanical models of crack propagation and their application to magma transport in dykes, J. Geophys. Res., 96, 10049–10077, 1991. 

Maccaferri, F., Bonafede, M., and Rivalta, E.: A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation, J. Volcanol. Geoth. Res., 208, 39–50, 2011. 

McKenzie, D.: The generation and compaction of partially molten rock, J. Petrol., 25, 713–765, 1984. 

McKenzie, D.: Constraints on melt generation and transport from U-series activity ratios, Chem. Geol., 162, 81–94, 2000. 

Minkowycz, W. J., Haji-Sheikh, A., and Vafai, K.: On departure from local thermal equilibrium in porous media due to a rapidly changing heat source: the Sparrow number, Int. J. Heat Mass Tran., 42, 3373–3385, 1999. 

Nield, D. A. and Bejan, A.: Convection in Porous Media, 3rd Edn., Springer Science+Business Media, Inc., ISBN 978-3-319-84189-2, 2006. 

Rivalta, E., Taisne, B., Buger, A. P., and Katz, R. F.: A review of mechanical models of dike propagation: Schools of thought, results and future directions, Tectonophysics, 638, 1–42, 2015. 

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

Rubin, A. M.: Propagation of magma-filled cracks, Annu. Rev. Earth Planet. Sc., 23, 287–336, 1995. 

Schmeling, H.: A model of episodic melt extraction for plumes, J. Geophys. Res., 111, B03202,, 2006. 

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. 

Schumann, T. E. W.: Heat transfer: A liquid flowing through a porous prism, J. Frankl. Inst., 208, 405–416, 1929. 

Spiegelmann, M., Kelemen, P. B., and Aharonov, E.: Causes and consequences of flow organization during melt transport: The reaction infiltration instability in compactible media, J. Geophys. Res., 106, 2061–2077, 2001. 

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 Tran., 24, 355–364,, 1981.  

Turcotte, D. and Schubert, G.: Geodynamics, Cambridge University Press, Cambridge, ISBN 978-0-521-18623-0, 2014. 

Verruijt, A.: Theory of Groundwater Flow, The Macmillan Press Ltd., London and Basingstoke, 141 pp.,, 1982. 

Wilcock, W. S. D. and Fisher, A. T.: Geophysical constraints on the subseafloor environment near Mid-Ocean ridges, 51–74, in: Subseafloor Biosphere, edited by: Cary, C., Delong, E., Kelley, D., and Wilcock, W. S. D., Washington DC, American Geophysical Union, ISBN 0-87590-409-2, 2004. 

Woods, A. W.: Flow in Porous Rocks: Energy and Environmental Applications, Cambridge University Press, Cambridge, 289 pp., ISBN 978-1-107-06585-7, 2015. 

Short summary
Fluid flow through rock occurs in many geological settings on different scales, at different temperature conditions and with different flow velocities. Fluid is either in local thermal equilibrium with the host rock or not. We explore the parameters of porous flow and give scaling laws. These allow us to decide whether porous flows are in thermal equilibrium or not. Applied to magmatic systems, moving melts in channels or dikes moderately to strongly deviate from thermal equilibrium.