Thermal nonequilibrium of porous flow in a resting matrix applicable to melt migration: a parametric study
 Institute of Geosciences, Goethe University, 60438 Frankfurt, Germany
 Institute of Geosciences, Goethe University, 60438 Frankfurt, Germany
Correspondence: Harro Schmeling (schmeling@geophysik.unifrankfurt.de)
Hide author detailsCorrespondence: Harro Schmeling (schmeling@geophysik.unifrankfurt.de)
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 nonequilibrium, the coupled heat equations for fluid and solid phases are formulated for a fluid migrating through a resting porous solid by porous flow. By nondimensionalizing the equations, two nondimensional numbers can be identified controlling thermal nonequilibrium: 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 nondimensional number, the initial thermal gradient G, which is the reciprocal of the nondimensional model height H. Three stages are observed: a transient stage followed by a stage with maximum nonequilibrium fluidtosolid temperature difference, ΔT_{max}, and a stage approaching the steady state. A simplified timeindependent ordinary differential equation for depthdependent (T_{f}−T_{s}) is derived and solved analytically. From these solutions simple scaling laws of the form $\left({T}_{\mathrm{f}}{T}_{\mathrm{s}}\right)=f\left(\mathit{Pe},\phantom{\rule{0.125em}{0ex}}G,z\right)$ are derived. Due to scaling they do not depend explicitly on ϕ anymore. The solutions for ΔT_{max} 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 ($>\mathrm{1}/G$) strong thermal nonequilibrium develops independently of Pe, (2) at low Pe ($<\mathrm{1}/G$) nonequilibrium decreases proportional to decreasing Pe⋅G, and (3) at low Pe (<1) and G of the order of 1 the scaling law is $\mathrm{\Delta}{T}_{max}\approx \mathit{Pe}$. The scaling laws are also given in dimensional form. The dimensional ΔT_{max} 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 nonequilibrium scale with the advective timescale in the highPe regime and with the interfacial diffusion time in the other two lowPe regimes. Applying the results to natural magmatic systems such as midocean ridges can be done by estimating appropriate orders of Pe and G. Plotting such typical ranges in the Pe–G 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 nonequilibrium with fluidtosolid temperature differences of up to several tens of kelvin, and (c) the dike regime is at full thermal nonequilibrium.
 Article
(8995 KB) 
Supplement
(486 KB)  BibTeX
 EndNote
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 supersolidus 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 subsolidus 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 nonequilibrium 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 viscoelastoplastic 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 subsolidus ambient rock, and thermal nonequilibrium delays freezing of the ascending melts and promotes initiation of further dikelike pathways. Indeed, for midoceanic ridges compositional nonequilibrium 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 nonequilibrium may also become important. Describing this nonequilibrium macroscopically, i.e., on a scale larger than the pores or channels, is the scope of this paper.
While the physics of thermal nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium of an arbitrary porous flow system depends on many parameters, our approach is to reduce the complexity of the system and systematically explore the nondimensional parameter space. It will be shown that only two nondimensional parameters control thermal nonequilibrium in porous flow, namely the Péclet number and the porosity. In our simple 1D model setup with constant flow velocity a third nondimensional number, the nondimensional initial thermal gradient G, is identified, which is equal to the reciprocal nondimensional model height $H=\mathrm{1}/G$. The nondimensionalization 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 nonequilibrium 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 midocean ridge setting.
2.1 Heat conservation equations
We start with considering a general twophase 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 nonequilibrium 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
For the definition of all quantities, see Table 1. Equation (1) can be rearranged into
Mass conservation for the fluid phase is given by
Inserting Eq. (3) into Eq. (2), conservation of energy for the fluid phase becomes
In a similar way, the conservation of energy of the solid phase is given by
which, assuming that v_{s}=0, is further simplified:
The term Q_{fs} 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:
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 Q_{fs} 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.
2.2 Scaling and nondimensionalization
Nondimensionalization is useful for interpreting models involving a large number of parameters. It usually helps reduce the number of parameters and identifies nondimensional parameters that control the evolution of the system. We write the two energy conservation equations in a nondimensional form, using
where ΔT_{0} 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; v_{f0} is the scaling fluid velocity; L is the scaling length
with ϕ_{0} as a scaling porosity; and t_{0} is the scaling time based on the diffusion time over the length L,
(see Table 1 for definitions). Primed quantities are nondimensional.
Introducing the fluidfilled pore width d_{f} and the solid width d_{s}, which may be the grain size or distance between fluid channels, the interfacial area density S scales with
for melt channels, tubes, pockets for all melt fractions, and melt films at small melt fractions, while S scales with
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 c_{s} 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
and
Equation (9a) shows that L scales both with the geometric mean of d_{f} and δ at small melt fractions and with the geometric mean of d_{s} and δ at large melt fractions. Thus, L is a natural length scale associated with thermal equilibrium of fluidfilled 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 d_{s} as the average distance between meltfilled pores or channels, which can be considerably larger than the grain size. Then both δ and d_{s} 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):
From Eqs. (4), (6), and (7) we get the nondimensional energy conservation equations for the fluid and solid phases, respectively:
From these equations we notice that the thermal evolution of the twophase system is controlled by two nondimensional numbers: the scaling porosity ϕ_{0} and the Péclet number Pe defined as
This number has already proven to be of high significance for determining whether thermal nonequilibrium 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 nondimensional, if not indicated otherwise.
In the following we consider a homogeneous twophase 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 nondimensional velocity v_{f}=1. This simplifies Eqs. (14) and (15) to
and
respectively. As we are interested in the evolution of the nonequilibrium temperature difference between the solid and fluid, subtraction of Eq. (18) from Eq. (17) gives
which is equivalent to
Note that while the temperatures T_{f} and T_{s} explicitly depend on two nondimensional numbers Pe and ϕ_{0}, the temporal evolution of the temperature difference (T_{f}−T_{s}) explicitly depends only on Pe. However, implicitly it is still a function of ϕ_{0} because T_{s} on the righthand side of Eq. (20) depends on ϕ_{0} via Eq. (18). Only for cases or stages with T_{s} independent of ϕ_{0} as proposed in Sect. 4 is the temperature difference (T_{f}−T_{s}) a function of only one nondimensional 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 nonequilibrium 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=\mathrm{1}/H$ is present in both phases (see Fig. 1). As a boundary condition both phase temperatures are set to 1 (nondimensional temperature) at z=0. At z=H a constant thermal gradient condition $\partial T/\partial z=\mathrm{1}/H$ (nondimensional) 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 T_{f}Peϕ_{0}, all nondimensional). This top boundary condition needs some justification: the hyperbolic partial differential Eqs. (17) or (18) require two welldefined 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.
This model setup adds a third nondimensional scaling parameter to the system, namely $G=\mathrm{1}/H$. It defines the initial nondimensional temperature gradient or conductive heat flux, positive for a flux directed upwards. To summarize, the temperatures depend on the nondimensional 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, $min\left(\mathrm{0.1},H/\mathrm{100}\right)$ for H<10. The time step was chosen as $\mathrm{d}t=\frac{\mathrm{1}}{\mathrm{4}}min\left(\mathrm{d}z/\mathit{Pe},\mathrm{d}{z}^{\mathrm{2}}\right)$, 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 ${\left(\right)}_{{\mathit{\delta}}_{q}=\frac{{q}_{\mathrm{tot}}\left(z=\mathrm{0}\right){q}_{\mathrm{tot}}\left(z=H\right)H\frac{\partial {T}_{\mathrm{mean}}}{\partial t}}{\left({q}_{\mathrm{tot}}\left(z=\mathrm{0}\right)+{q}_{\mathrm{tot}}\left(z=H\right)\right)/\mathrm{2}}}max$ , where q_{tot} is the total nondimensional vertical heat flux (conductive and advective) and T_{mean} 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.
First, some example numerical results are shown in Fig. 2 to understand the physics and the typical behavior.
3.1 Evolution of temperatures and thermal nonequilibrium 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 T_{f} and T_{s} 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 depthdependent conductive and advective heat fluxes through the solid and fluid phases, respectively. Figure 2e shows the evolution of T_{f} and T_{s} 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 (T_{f}−T_{s}) 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 nonzero; 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 steadystate 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 T_{s} 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 T_{s} is no longer constant in time, different melt fractions lead to different rates of temperature increase (Fig. 2e) and also to different evolutions of (T_{f}–T_{s}) (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 T_{f} 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 (T_{f}−T_{s}).
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 firstorder estimate of thermal nonequilibrium 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 nondimensional parameters Pe and $G=\mathrm{1}/H$, and ΔT_{max} has been determined for each model (Fig. 3). Some first observations can be made.

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

ΔT_{max} is proportional to G, i.e., to the nondimensional temperature gradient for G<0.1.

ΔT_{max} 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.

ΔT_{max} 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 ΔT_{max} could be derived, based on the two scaling parameters. In the next section, we propose an analytical derivation of ΔT_{max} values to obtain scaling laws and confirm the observed proportionalities.
In this section a simplified analytical solution for the zdependent temperature difference between fluid and solid will be derived. From this solution the maximum temperature differences ΔT_{max} 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 nonequilibrium temperature difference (T_{f}−T_{s}). We simplify the problem by considering the hypothetical case in which (T_{f}−T_{s}) does not change with time and, moreover, in which the thermal gradient in the solid phase is fixed and linear, with $\partial {T}_{\mathrm{s}}/\partial z=G=\mathrm{1}/H$ (nondimensional, with dimensions $G=\mathrm{\Delta}{T}_{\mathrm{0}}/H)$. Although different from initial and steadystate 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 T_{f} and T_{s} was observed to be quite similar indeed. Besides, at the end of stage 1 (Sect. 3.1), T_{s} remains close to initial conditions; therefore a fixed linear gradient of slope $G=\mathrm{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
While in the general case of Eq. (20) the temperature difference implicitly depends on ϕ_{0}, i.e., on the three nondimensional parameters Pe, ϕ_{0}, and G, Eq. (21) no longer depends on ϕ_{0} because we replaced $\partial {T}_{\mathrm{s}}\left({\mathit{\varphi}}_{\mathrm{0}}\right)/\partial z$ with −G, which is independent of ϕ_{0}. Equation (21) is a secondorder ordinary differential equation for (T_{f}−T_{s}) whose solution can be analytically derived as (see Supplement for details)
where r_{1} and r_{2} are the roots of the associated equation of Eq. (21):
The parameters α and β are constrained by the boundary conditions $\left({T}_{\mathrm{f}}{T}_{\mathrm{s}}\right)=\mathrm{0}$ at z= 0 and $\frac{\partial \left({T}_{\mathrm{f}}{T}_{\mathrm{s}}\right)}{\partial z}=\mathrm{0}$ at z=H:
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 depthdependent temperature difference (T_{f}−T_{s}) can be determined. It is found that the maximum is always at z=H. This value will be denoted as ΔT_{max} 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 ΔT_{max} 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 timeindependent Eq. (21) to obtain an analytical solution of an intrinsically timedependent process as long as we are interested only in the maximum value of (T_{f}−T_{s}). 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 $\mathrm{\Delta}{T}_{max}=\mathrm{1}$ where the ΔT_{max}(Pe) curves become nonlinear (Fig. 3a).
4.3 Scaling laws for temperature differences at certain parameter limits
The analytical solution for ΔT_{max} 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 nonequilibrium 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 (T_{f}−T_{s}), i.e., also ΔT_{max}, 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
With this condition Eq. (22) tends to the following limit (see Supplement):
with
which simplifies for $z=H=\mathrm{1}/G$ to
This is the limit for Pe→0. This limit gives predictions for ΔT_{max} 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 $\mathit{Pe}<\mathrm{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
Equation (29) confirms the proportionalities observed in Fig. 3, namely $\mathrm{\Delta}{T}_{max}\propto \mathit{Pe}$ (Fig. 3a) and $\mathrm{\Delta}{T}_{max}\propto G$ (Fig. 1b).
4.3.2 Limit Pe→∞
To obtain the limit of Eq. (22) for Pe→∞, Eq. (22) can be linearized with respect to $\mathrm{4}/{\mathit{Pe}}^{\mathrm{2}}\ll \mathrm{1}$. Applying the rule of L'Hôpital, Eq. (22) tends to the following limit:
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 ΔT_{max} 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 midocean 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 10^{5} in dike systems (Schmeling et al., 2018).
While the melt fraction does not influence ΔT_{max} (cf. Eqs. 22, 30), it does influence the longterm temporal behavior because T_{s} 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=\mathrm{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 d_{f} and the interfacial boundary layer thickness δ (Eqs. 9 with 11). L would evolve nonlinearly 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 nondimensional G ranges between an order of 1 and an order of 10^{−5}.
In Fig. 4 we explore ΔT_{max} variations using the analytical solution Eq. (22), in which ΔT_{max} depends on Pe and G. Three main regimes can be distinguished.

Regime 1. For high Pe values, (T_{f}−T_{s}) tends to the relationship described in Eq. (30). The temperature difference increases linearly with distance from the bottom (z=0), reaching $\mathrm{\Delta}{T}_{max}=\mathrm{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 ΔT_{max} to G disappears because the maximum value of z is equal to $H=\mathrm{1}/G$.

Regime 2. For Pe≪1, or more precisely, for $\mathit{Pe}\ll \frac{\mathrm{1}}{G}$ represented by the oblique dashed line in Fig. 4, (T_{f}−T_{s}) varies with distance from the bottom according to $\left(\mathrm{1}{e}^{z}\right)$ and is proportional to Pe and G. This means that large temperature gradients favor large temperature differences. In this domain, (T_{f}−T_{s}) 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, (T_{f}−T_{s}) tends to the relationship proposed in Eq. (26). In this domain, (T_{f}−T_{s}) 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 nonlinearly from about 0 to 0.4 with increasing z.
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 ΔT_{max} 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 timedependent temperature differences (T_{f}−T_{s}) 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 longterm behavior is approached. We emphasize that this analytical solution is a very good approximation of the depthdependent temporal maximum temperature difference that can be reached in such porous systems.
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 (T_{f}−T_{s}) 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 onesided (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 buildup of T_{f}−T_{s} 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 $\left(\mathit{Pe}{T}_{\mathrm{influx}}\mathit{Pe}{T}_{\mathrm{outflux}}=\mathit{Pe}\right)$ into a high conductive outflux $\left(\partial T/\partial z\right)$ at the top. It is unlikely that such situations occur in natural systems.
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 nondimensional parameter, the lid thickness H_{lid}. 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 (T_{f}−T_{s}) values occur near the top, but they disappear for H_{lid}≥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 (T_{f}−T_{s}) values than in Fig. 6 (top row) close to the top. However, typical natural magmatic systems are on the H_{lid}>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 T_{f}−T_{s} 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 nonlinear initial temperature profile between ${T}_{\mathrm{f}}={T}_{\mathrm{s}}=\mathrm{1}$ at the bottom and ${T}_{\mathrm{f}}={T}_{\mathrm{s}}=\mathrm{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 nonequilibrium (see e.g., Eq. 22, which explicitly contains the temperature gradient G), the above results are expected to be different, and a stronger thermal nonequilibrium is expected in regions with higher gradients. Schmeling et al. (2018) used a step function with ${T}_{\mathrm{f}}={T}_{\mathrm{s}}=\mathrm{1}$ at z=0 and ${T}_{\mathrm{f}}={T}_{\mathrm{s}}=\mathrm{0}$ at z>0 as the initial condition, i.e., an extremely nonlinear 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 timeindependent 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 timedependent 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 nonequilibrium compared to settings with nonlinear initial temperature profiles.
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 nondimensional numbers are introduced, we nondimensionalize 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 $\stackrel{\mathrm{\u0303}}{L}=\frac{L}{\sqrt{{\mathit{\lambda}}_{\mathrm{eff}}^{\prime}}}$ with ${\mathit{\lambda}}_{\mathrm{eff}}^{\prime}={\mathit{\lambda}}_{\mathrm{eff}}/{\mathit{\lambda}}_{\mathrm{f}}$ and $\stackrel{\mathrm{\u0303}}{\mathit{Pe}}=\frac{\mathit{Pe}}{\sqrt{{\mathit{\lambda}}_{\mathrm{eff}}^{\prime}}}$. With this scaling, Eqs. (14) and (15) turn into (for clarity, primes indicate nondimensional quantities)
and
Inspection of these equations shows that two more nondimensional numbers are introduced: the ratio of diffusivities ${\mathit{\kappa}}_{\mathrm{s}}^{\prime}$ and the ratio of the products' density and heat capacity, ${\mathit{\rho}}_{\mathrm{s}}^{\prime}{c}_{p,\mathrm{s}}^{\prime}$.
As Eqs. (31) and (32) cannot be merged into one timeindependent ordinary differential equation for (T_{f}−T_{s}) as in Sect. 4.1, we numerically tested some cases with $\stackrel{\mathrm{\u0303}}{\mathit{Pe}}=\mathit{Pe}=\mathrm{1}$ and ${\mathit{\lambda}}_{\mathrm{eff}}^{\prime}=\mathrm{1}$ in which the diffusivity ratio and the ratio of ${\mathit{\rho}}_{\mathrm{s}}^{\prime}{c}_{p,\mathrm{s}}^{\prime}$ were varied between 0.1 and 10 (see Fig. 8). The results show that for the fixed combination of Pe=1 and ${\mathit{\lambda}}_{\mathrm{eff}}^{\prime}=\mathrm{1}$ the magnitude of thermal nonequilibrium 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 ${\mathit{\kappa}}_{\mathrm{s}}^{\prime}=\mathrm{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 ${\mathit{\kappa}}_{\mathrm{s}}^{\prime}=\mathrm{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., ${\mathit{\rho}}_{\mathrm{s}}^{\prime}{c}_{p,\mathrm{s}}^{\prime}$, Fig. 8e and f show that a high value slows down the longterm timedependent variations, while a small value leads to rapid longterm temporal variations in (T_{f}−T_{s}) and faster convergence to the steady state, which is similar to the case with equal properties.
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 ${\mathit{\rho}}_{\mathrm{s}}^{\prime}{c}_{p,\mathrm{s}}^{\prime}$ 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., ${\mathit{\kappa}}_{\mathrm{s}}^{\prime}$ 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 ${\mathit{\lambda}}_{\mathrm{eff}}^{\prime}=\mathrm{1}$) equal to the cases depicted in Fig. 5. The timedependent 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 nonequilibrium temperature differences shown in the regime diagram Fig. 4 are also applicable to a water–sedimentary rock system.
5.2 Timescales
It is interesting to evaluate the timescales for reaching the maximum nonequilibrium 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, t_{90 %}, and the time needed to reach steady state, t_{steady}. The latter has been determined as the time at which the maximum difference between (T_{f}(z)−T_{s}(z)) curves at two subsequent time steps becomes less than ${\mathrm{10}}^{\mathrm{8}}\mathrm{\Delta}{T}_{max}$. 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 d_{char}, giving ${t}_{\mathrm{advd}}={d}_{\mathrm{char}}/{v}_{\mathrm{f}\mathrm{0}}$, or on diffusion over the characteristic distance, giving ${t}_{\mathrm{diffd}}={d}_{\mathrm{char}}^{\mathrm{2}}/{\mathit{\kappa}}_{\mathrm{0}}$. We tested these timescales with the two natural length scales of the models. The first is the scaling length L (equal to 1 nondimensional), representing essentially the geometric mean of the channel width of the pores, d_{f}, 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 t_{90 %} and t_{steady} 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).
In regime 1 (high Pe), t_{90 %} is proportional to t_{advH} (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 t_{steady}∝t_{diffH}. 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 nonequilibrium is still present.
In regime 2 (low Pe and G<0.1, i.e., H>10) the time for reaching ΔT_{max} is controlled by interfacial heat transfer (Fig. 10a, red asterisks) on the length scale L, resulting in t_{90 %} proportional to t_{0}. 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 ΔT_{max} 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 nonequilibrium 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 t_{steady}∝t_{diffH}. Clearly, it is also controlled by diffusion.
5.3 Applications to magmatic systems
We now test the possible occurrence of thermal nonequilibrium in natural magmatic systems based on the suggested controlling nondimensional parameters, namely the Péclet number Pe, the initial thermal gradient G ($=\mathrm{1}/H$), and the melt fraction ϕ. Typical stages of melt flow for midocean 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 highporosity (> 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 Darcyflowrelated Péclet number:
As we preferably use the melt pore dimension d_{f} in our scalings (Eqs. 9a and 10a), we need to relate it to the solid phase dimension d_{s} by using
Using Eqs. (34), (9a), and (16), we arrive at the Péclet number used here.
Schmeling et al. (2018) reviewed and estimated typical pore or channel spacings d_{s} 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 Pe_{D} 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.
>10^{5} 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 timedependent, a good estimate is $\mathit{\delta}={c}_{\mathrm{th}}\sqrt{{\mathit{\kappa}}_{\mathrm{0}}t}$ (in dimensional form), where c_{th} 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 v_{0} and system height H, i.e., by $t=H/{v}_{\mathrm{0}}=H\mathit{\varphi}/{v}_{\mathrm{D}}$, we may express v_{D} in terms of the Péclet number Pe_{D}. With the resulting t and δ we arrive at the following Péclet number (H and d_{f} are dimensional or nondimensional):
For midocean 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 d_{f} 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 10^{4} to 10^{7} for the dike stage (c). To estimate typical nondimensional thermal gradients G^{′} (or layer thickness H^{′}), the above estimate for δ and d_{f} can be inserted into the scaling length L (Eq. 9a) to arrive at a nondimensional ${G}^{\prime}=\mathrm{1}/{H}^{\prime}$.
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 nonequilibrium 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 nonequilibrium at sufficiently high thermal gradients.
 c.
After transition to diking full thermal nonequilibrium 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 midocean ridges. Thus, Péclet numbers are expected to be smaller, but nondimensional 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 nonequilibrium 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
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 d_{f} or d_{s}, 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 nonequilibrium temperature difference is simply equal to the imposed temperature difference, $\mathrm{\Delta}{T}_{max}=\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{T}_{\mathrm{0}}$, and the characteristic time to reach maximum nonequilibrium is simply ${t}_{\mathrm{char}}=H/{v}_{\mathrm{f}\mathrm{0}}$, i.e., the total time of a fluid particle for passing through the system.

Regimes 2 and 3. For small Péclet numbers ($\mathit{Pe}<\frac{H\sqrt{S}}{\sqrt{{\mathit{\varphi}}_{\mathrm{0}}\left(\mathrm{1}{\mathit{\varphi}}_{\mathrm{0}}\right)\mathit{\delta}}}$) the maximum temperature difference scales like
$$\begin{array}{}\text{(39)}& \begin{array}{rl}\mathrm{\Delta}{T}_{max}& ={\displaystyle \frac{G{v}_{\mathrm{f}\mathrm{0}}{\mathit{\varphi}}_{\mathrm{0}}\left(\mathrm{1}{\mathit{\varphi}}_{\mathrm{0}}\right)\mathit{\delta}}{{\mathit{\kappa}}_{\mathrm{0}}S}}\\ & ={\displaystyle \frac{{v}_{\mathrm{f}\mathrm{0}}{\mathit{\varphi}}_{\mathrm{0}}\left(\mathrm{1}{\mathit{\varphi}}_{\mathrm{0}}\right)\mathit{\delta}}{H{\mathit{\kappa}}_{\mathrm{0}}S}}\mathrm{\Delta}{T}_{\mathrm{0}},\end{array}\end{array}$$and the characteristic time for reaching this nonequilibrium scales with t_{0}, i.e.,
$$\begin{array}{}\text{(40)}& {t}_{\mathrm{char}}={\displaystyle \frac{{\mathit{\varphi}}_{\mathrm{0}}\left(\mathrm{1}{\mathit{\varphi}}_{\mathrm{0}}\right)\mathit{\delta}}{{\mathit{\kappa}}_{\mathrm{0}}S}}.\end{array}$$These relations can easily be used to assess the potential of thermal nonequilibrium in systems of fluid flow through solids with given geometrical properties and fluid fractions.
In the above discussion we used the terms moderate thermal nonequilibrium, which we may identify with $\mathrm{\Delta}{T}_{max}^{\prime}$ of a few percent to, say, 30 %, while full thermal nonequilibrium includes higher values up to 100 %. To translate this into dimensional ΔT_{max} values, what are typical ΔT_{0} values for midocean 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 T_{asth} and, say, 0.8T_{asth}, i.e., $\mathrm{\Delta}{T}_{\mathrm{0}}=\mathrm{0.2}{T}_{\mathrm{asth}}\cong \mathrm{200}$ K. Thus moderate nonequilibrium may be represented by excess temperatures of the melt with respect to the solid between, say, 6 and 60 K, while full thermal nonequilibrium suggests the full ΔT_{0} range of 60 to 200 K or even higher in the case of dikes extending up into the lithosphere with temperatures below 0.8 T_{asth}. 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.
In conclusion, we showed that in magmatic systems characterized by twophase flows of melts with respect to solid, thermal nonequilibrium between melt and solid may arise and becomes important under certain conditions. The main conclusions are summarized as follows.
From nondimensionalization of the governing equations, three nondimensional numbers can be identified controlling thermal nonequilibrium: the Péclet number Pe, the melt porosity ϕ, and the initial nondimensional temperature gradient G in the system. The maximum possible nonequilibrium solid–fluid temperature difference ΔT_{max} is controlled only by two nondimensional numbers: Pe and G. Both numerical and analytical solutions show that in a Pe–G parameter space three regimes can be identified.

In regime 1 (high $\mathit{Pe}(>\mathrm{1}/G)$) strong thermal nonequilibrium develops independently of Pe, and a nondimensional scaling law ${T}_{\mathrm{f}}{T}_{\mathrm{s}}=Gz$ has been derived.

In regime 2 (low $\mathit{Pe}(<\mathrm{1}/G)$ and low G(<0.3)) nonequilibrium decreases proportionally to decreasing Pe and G, and the nondimensional scaling law reads ${T}_{\mathrm{f}}{T}_{\mathrm{s}}=\mathit{Pe}G\left(\mathrm{1}{e}^{z}\right)$.

In regime 3 (low Pe(<1) and G of the order of 1) nonequilibrium scales with Pe and G and is depthdependent. The scaling law is ${T}_{\mathrm{f}}{T}_{\mathrm{s}}=\mathit{Pe}G\phantom{\rule{0.125em}{0ex}}\left(\mathrm{1}M\left(z\right)\right)$, where M(z) depends on G.
Further conclusions include the following.

The timescales for reaching thermal nonequilibrium scale with the advective timescale in the highPe regime and with the interfacial diffusion time in the other two lowPenumber regimes.

Applying the results to natural magmatic systems such as midocean ridges can be done by estimating appropriate orders of Pe and G. Plotting such typical ranges in the Pe–G 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 nonequilibrium, and (c) the dike regime is at full thermal nonequilibrium.

In the studied setup G was constant, leading to conservative estimates of thermal nonequilibrium. Any other depthdependent initial temperature distributions generate higher nonequilibrium than reported here.

The derived scaling laws for thermal nonequilibrium are valid for equal solid and fluid properties. Assuming different properties such as for a water–sandstone system results in similar maximum nonequilibrium 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 nonequilibrium for anastomosing systems (Hart, 1993).
The MATLAB code is listed in the Supplement and is available on request.
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: https://doi.org/10.5194/se1310452022supplement.
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.
The contact author has declared that neither they nor their coauthor 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.
This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. 403710316).
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 NonThermal Equilibrium, NonDarcian, Variable Porosity Incompressible 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 nonequilibrium in heterogeneous media, Springer Science+Business Media, Inc., https://doi.org/10.1007/9783319146669, 2016.
Furbish, D. J.: Fluid Physics in geology, Oxford University Press, New York, 476 pp., ISBN 0195077016, 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, https://doi.org/10.1029/94JB02544, 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, https://doi.org/10.1016/00179310(94)903581, 1994.
Landwehr, D., Blundy, J., ChamorroPerez, E. M., Hill, E., and Wood, B.: Useries disequilibria generated by partial melting of spinel lherzolite, Earth Planet. Sc. Lett., 188, 329–348, 2001.
Lister, J. R. and Kerr, R. C.: Fluidmechanical 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 Useries activity ratios, Chem. Geol., 162, 81–94, 2000.
Minkowycz, W. J., HajiSheikh, 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 9783319841892, 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 melttransport: Implications for the evolution of the lithosphereasthenosphere boundary, arXiv [preprint], arXiv:2009.01496, 2020.
Rubin, A. M.: Propagation of magmafilled 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, https://doi.org/10.1029/2004JB003423, 2006.
Schmeling, H., Marquart, G., and Grebe, M.: A porous flow approach to model thermal nonequilibrium 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, https://doi.org/10.1016/00179310(81)900430, 1981.
Turcotte, D. and Schubert, G.: Geodynamics, Cambridge University Press, Cambridge, ISBN 9780521186230, 2014.
Verruijt, A.: Theory of Groundwater Flow, The Macmillan Press Ltd., London and Basingstoke, 141 pp., https://doi.org/10.1007/9781349167692, 1982.
Wilcock, W. S. D. and Fisher, A. T.: Geophysical constraints on the subseafloor environment near MidOcean 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 0875904092, 2004.
Woods, A. W.: Flow in Porous Rocks: Energy and Environmental Applications, Cambridge University Press, Cambridge, 289 pp., ISBN 9781107065857, 2015.
 Abstract
 Introduction
 Governing equations and model setup
 Numerical model results
 Scaling laws derived from analytical solution
 Discussion
 Conclusions
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Governing equations and model setup
 Numerical model results
 Scaling laws derived from analytical solution
 Discussion
 Conclusions
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement