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

. 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 Peclet number 𝑃𝑒 10 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, ∆𝑇 𝑚𝑎𝑥 , and a stage approaching the steady state. A simplified time-independent ordinary differential equation for depth-dependent (𝑇 𝑓 − 𝑇 𝑠 ) is 15 derived and solved analytically. From these solutions simple scaling laws of the form (𝑇 𝑓 − 𝑇 𝑠 ) = 𝑓(𝑃𝑒, 𝐺, 𝑧) are derived. Due to scaling they don’t depend explicitly on 𝜙 anymore. The solutions for ∆𝑇 𝑚𝑎𝑥 and the scaling laws are in good agreement with the numerical solutions. The parameter space 𝑃𝑒, 𝐺 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 𝑃𝑒 ∙ 𝐺 ; 3) at low Pe (<1) and G of order 1 the scaling law is ∆𝑇 𝑚𝑎𝑥 ≈ 𝑃𝑒 . The scaling laws are also 20 given in dimensional form. The dimensional ∆𝑇 𝑚𝑎𝑥 depends on the initial temperature gradient, the flow velocity, the


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 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 also water may 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 dykes 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 melts moving in dykes through a cold elastic or visco-elasto-plastic ambient rock, the fluids are assumed as isothermal (e.g.Maccaferri et al., 2011;Keller et al., 2013).However, on local scale of channel or dyke 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 dyke 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 dyke-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.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 general, on the Peclet 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 nondimensional parameters control thermal non-equilibrium in porous flow, namely the Peclet number and the porosity.In our simple 1D model setup with constant flow velocity a third non-dimensional number, the model height H=1/G, where G is the non-dimensional initial thermal gradient is identified.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 to an anastomosing melt ascent system typical for mid-oceanic ridges in a second paper (Chevalier and Schmeling, in prep).

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 e.g. by de Lemos (2016).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 (3) into (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   = 0, is further simplified: The term   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, 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   and allows for a 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.

Scaling and non-dimensionalization
Non-dimensionalization is useful for interpreting models involving a large number of parameters.It usually helps reducing 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 where  0 is the macroscopic scaling temperature difference of the system, i.e. the initial temperature difference between top and bottom, x,y,z is a distance, vf0 is the scaling fluid velocity, L is the scaling length with  0 as a scaling porosity, and  0 is the scaling time based on the diffusion time over the length , (see Table 1 for definitions).Primed quantities are non-dimensional.Introducing the fluid filled pore width   and the solid (grain) width ds, the interfacial area density S scales with for melt channels, tubes, pockets for all melt fractions, and for 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 2 for melt channels, 4 for melt tubes, 6 for melt pockets, and 2 for melt films at small melt fractions.The geometrical constant cs is of order 2 for melt channels, and 6 for melt films or suspensions.Thus, the scaling time and scaling length can also be written as and Eq. (9a) shows that L scales both with the geometric mean of   and  at small melt fractions, and with the geometric mean of   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.
We assume that the fluid and solid phases have the same densities and thermal properties (but relax this assumption later in section 5.1.3): , =  , =  ,0 ,   =   =  0 ,   =   =    ,0  0 =  0 (13) From Eq. ( 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 Peclet number Pe defined as 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   = 1.This simplifies equations ( 14) and (15) to and 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: 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 (  −   ) 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 section 4, the temperature difference (  −   ) is a function of only one non-dimensional parameter, Pe, and no more of  0 .

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 −  = −1  ⁄ is present in both phases (see Fig. 1).As boundary condition both phases temperatures are set to 1 (non-dimensional temperature difference) at z = 0.At z = H a constant thermal gradient condition    ⁄  = − 1  ⁄ (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   0 , while at z = H it evolves freely with the fluid temperature, i.e. it is given by     0 (all non-dimensional).This top boundary condition needs some justification: The hyperbolic partial differential equations Eq. ( 17 This model setup adds a third non-dimensional scaling parameter to the system, namely  = 1/.It defines the initial nondimensional 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.

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  = 0.1 or (0.1,/100) for  < 10.The the time step was chosen as  = 1 4 (/, 2 ), i.e. taking the minimum of the Courant or diffusion criterion.
Tests with smaller spatial and temporal resolution have been carried out and did not change the results visibly.

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

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: H = 10, ϕ = 0.1, model 2: forcing the solid temperature to progressively increase.At the end of stage 1 the maximum temperature difference is approached (Fig. 2f).Because the solid temperature hasn't risen significantly at that time (at t = 4 in the example) compared to the fluid temperature (Fig. 2e) different melt fractions do not affect the temperature differences during this stage (Fig. 2f in which all curves merge in 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. 2c), the temperature difference remains constant and at maximum at the top (Fig. 2d).Solid-fluid heat transfer is at maximum during this stage.As Ts is no more constant in time, different melt fractions lead to different rates of temperature increase (Fig. 2c) and also to different evolutions of (Tf -Ts) (Fig. 2d solid curves compared to dashed curves).A higher melt fraction increases the heat transfer into the solid (c.f. last term in Eq. 18), resulting in a faster increase of the solid temperature whose gradient flattens earlier.Thus, the end of stage 2 is reached earlier (Fig. 2b).
Stage 3: As the fluid temperature rises close to the Tf value at the bottom, its increase slows down, and heat transfer, 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 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 of 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 (stage1), here after t = 4 (Fig. 2f).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. 2d).The higher in the model, the longer the temperature difference remains at maximum.A higher melt fraction accelerates the decrease of (Tf -Ts).The absolute maximum temperature difference in space and time does not depend on boundary conditions (see also section 5.1.2where the influence of boundary conditions is discussed), nor on the z-position nor 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 nonequilibrium conditions and possible temperature difference in a magmatic system.In the following sections we study how this maximum temperature difference evolves when varying the parameter Pe.

Maximum temperature difference
The maximum temperature difference of a model can be defined as the maximum value reached in space and time (c.f.Fig. 2d).A series of models has been carried out for the two different non-dimensional parameters Pe, and  = 1/, and ΔTmax has been determined for each model (Fig. 3).Some first observations can be made:  For all Pe, Δ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 order 1, i.e. when H reaches 1 or the dimensional H reaches the scale L.
 ΔTmax is essentially independent of ϕ as models with different ϕ 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.

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.

Analytical solution of the governing equations
We are interested in an analytical solution of the equation ( 20) controlling the non-equilibrium temperature difference (  −   ).We simplify the problem by considering the hypothetical case in which (  −   ) does not change with time, and, moreover, in which the thermal gradient in the solid phase is fixed and linear, with     ⁄  = − = − 1  ⁄ (nondimensional, with dimensions:  =  0 /).Although different from initial and steady-state stages described in the 1D models (section 3.1), this hypothetical case is quite similar to what can be observed at the very beginning of the second stage described in section 3.1 (c.f.Fig. 2d,f).In this second stage, the evolution of Tf and Ts was observed being quite similar indeed.Besides, at the end of stage 1 (section 3.1), Ts remains close to initial conditions, therefore a fixed linear gradient of slope − = −1/ is justified.Since the maximum temperature difference between the two phases is observed starting from the end of stage 1 and during stage 2 (section 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 non-dimensional parameters Pe,  0 , and G, Eq. ( 21) does no more depend on  0 because we replaced    ( 0 )  ⁄  by -G which is independent of  0 .Eq. ( 21) is a second order ordinary differential equation for (  −   ) whose solution can be analytically derived as (see supplementary material for details) where r1 and r2 are the roots of the associated equation of Eq. ( 21) The parameters α and β are constrained by the boundary conditions: (  −   ) = 0 at z = 0 and The third term in Eq. ( 22) is a particular solution for Eq. ( 21).

Comparison with numerical models
From Eq. ( 22) the maximum value of the depth-dependent temperature difference (  −   ) 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 the differences become > 1 %, up to 6%.This general good agreement is another justification for using the time-independent equation ( 21) to obtain an analytical solution of an intrinsically time-dependent process as long as we are interested only in the maximum value of (  −   ).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).

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 on the relative influences of the two controlling parameters  and G, described in section 2.2 and 2.3.The Peclet number is already known to be of great importance for thermal equilibrium/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.
Eq. ( 22) is, however, complicated, and the assessment of the relative importance of  and G for different possible regimes is limited.In this section, we study the evolution of (  −   ), i.e. also ΔTmax, in a few limiting cases.This enables us better understanding each parameter influence and to derive some scaling laws for different regimes.

Limit 𝑷𝒆 → 𝟎
When Pe tends to 0, we have the condition With this condition Eq. ( 22) tends to the following limit (see supplementary material): with which simplifies for  =  = 1/ to This is the limit for Pe  0. This limit gives predictions for ΔTmax in very good agreement with Eq. ( 22) for Pe < 1 (having  = 0.1) (see Fig. S1 in the supplementary material).In the limit  → 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 Eq. ( 29) confirms the proportionalities observed in Fig. 3, namely ∆  ∝  (Fig. 3a), and ∆  ∝  (Fig. 1b), respectively.

Limit Pe  ∞
To obtain the limit of Eq. ( 22) for  → ∞, Eq. ( 22) can be linearized with respect to 4  2 ⁄ ≪ 1. Applying the rule of L'Hôpital Eq. ( 22) tends to the following limit: For details, see supplementary material.This limit is also the solution of Eq. ( 21) when neglecting the diffusive and heat transfer terms.As demonstrated in the supplementary material this limit predicts ΔTmax values in very good agreement with Eq. ( 22) for Pe > 100.

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 order 10 -5 to 10 -3 in partially molten regions with distributed porous flow to higher values of order 1 or larger at depths where channels have merged, and further to very high values of order 10 5 in dyke systems (Schmeling et al., 2018).
While the melt fraction does not influence ∆  (c.f.Eq. ( 22, 30)) it influences the long term temporal behavior once Ts is  0dependent (c.f.Eq. ( 20).Therefore 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, nor in dyke 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 dimension  = 1/ 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  (Eq. 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 dykes.As will be shown in section 5.3 in more detail the resulting non-dimensional G ranges between order 1 to order 10 -5 .
In Figure 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, (  −   ) 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 maximum 1 while the solid temperature increases linearly with z from 0 to 1.
 Regime 2: For  ≪ 1, or more precisely, for  ≪ 1  represented by the oblique dashed line in Fig. 4, (  −   ) varies with distance from the bottom according to (1 −  − ), and is proportional to Pe and G.This means that large temperature gradients favor large temperature differences.In this domain, (  −   ) tends to the relationships presented in Eq. ( 29).
 Regime 3: For large initial temperature gradient G close to 1 (small H) and  ≪ 1, (  −   ) tends to the relationship proposed in Eq. ( 26).In this domain, (  −   ) is proportional to Pe but no more to G. The depth-dependence is given by (1 − ()) which at G = 1 increases non-linearly from about 0 to 0.4 with increasing z.

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 very well the numerical results.This is shown in Fig. 5 where for various combinations of Pe and G the time-dependent temperature differences (  −   ) 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 the 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.

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 (regime 2 and 3, low Pe), and for early stages for regime 1 but might be inappropriate for high temperature increases (high Peregime 1) at later stages (see section 4.3.4).In order to quantify the influence of this choice of boundary conditions on our results, we compared the evolution of (  −   ) -profiles for three Peclet numbers and two heights H, 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.
 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 equations ( 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 section 4.1.Numerically it works well for our system without producing instabilities or oscillations.
Comparing the top and bottom row of Fig. 6, the constant temperature gradient condition produces quite similar results as the open boundary condition for all Pe and H values tested during the first and second stage of temporal evolution (c.f.section 3.1).The agreement becomes worse for stage 3 when approaching steady state at large Pe.Comparing the other two boundary conditions (2 nd and 3 rd row of Fig. 6) with the constant gradient condition (top row) shows that the effect of the top boundary during stage 1 and 2 is still small sufficiently far away from the top.Only for the small Pe -case (left column of Fig. 6) 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 high H with zero temperature boundary condition (3 rd row 4 th column in Fig. 6) shows a strong build-up of   −   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 (   −    = ) into a high conductive outflux (/) at the top.It is unlikely that such situations occur in natural systems.
In summary, the influence of boundary conditions on fluid and solid temperatures evolution depends mostly on the domain size (H) and on the value of Pe.The larger these two parameters, the less important is the influence of boundary conditions within almost the whole model domain.If one is interested in the maximum value of   −   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   =   = 1 at the bottom and   =   = 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 ), 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   =   = 1 at z = 0 and   =   = 0 at z > 0 as initial condition, i.e. an extremely non-linear profile near z = 0. Assuming this initial temperature profile Figure 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 case (Eq.22) is also shown.As expected, at early stages the temperature differences are significantly larger than given by the analytical solutions by a factor 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 non-equilibrium compared to settings with non-linear initial temperature profiles.

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.Eq. ( 14) and ( 15) turn into (for clarity, primes indicate non-dimensional quantites): and Inspection of these equations shows that three more non-dimensional numbers are introduced: the ratio of diffusivities   ′, the ratio of the products density and heat capacity,   ′ , ′, and a new effective conductivity for heat transfer,   ′.
As equations ( 32) and ( 33) cannot be merged into one time-independent ordinary differential equation for (  −   ) as in section 4.1, we numerically tested some cases with  = 1 and   ′ = 1 in which the diffusivity ratios and the ratios of   ′ , ′ were varied between 0.1 and 10 (see Fig. 8).The results show that for the fixed combination of  = 1 and   ′ = 1 the magnitude of thermal non-equilibrium remains in the same order of magnitude (0.1) as for equal properties (Fig. 8).The time-dependence is significantly affected: For a high ratio of   ′ = 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   ′ = 0.1 the solid temperature departs more strongly from the initial linear gradient, and the solidfluid temperature difference slowly drops with time on the long term.
Varying the potential to store heat in the solid, i.e.   ′ , ′, Fig. 8e and f shows that a high value slows down the long term time-dependent variations, while a small value leads to rapid long term temporal variations of (  −   ) and faster convergence to the steady state which is similar to the equal properties case.
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 three times larger than that of rock, and the density is almost three times less, the product   ′ , ′ is about 0.78, i.e. of order 1.However, the thermal diffusivity of water is significantly smaller than that of rock, typically by a factor 16, i.e.   ′ is about 16.We tested a few cases (Fig. 9) with Peclet numbers and initial thermal gradients  (i.e.inverse model heights) (assuming for simplicity   ′ = 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.

Time scales
It is interesting to evaluate the time scales 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,  90% , and the time needed to reach steady-state,   .The latter has been determined as the time at which the maximum difference between (  () −   ())curves at two subsequent time steps becomes less than 10 −8 ∆  .These times can be compared with different time scales that may characterize the evolution of temperatures in the models.These time scales can be based on advection over a characteristic distance  ℎ giving   =  ℎ  0 ⁄ , or on diffusion over the characteristic distance giving   =  ℎ 2  0 ⁄ .We tested these time scales with the two natural length scales of the models.
The first is the scaling length L (= 1 non-dimensional) representing essentially the geometric mean of the channel width of the pores,   , and the interfacial boundary layer thickness .The second is the model height .Grouping the models depending on the regime they belong to (see section 4.3.4,and Fig. 4), we plotted the recorded times  90% and   versus the characteristic time scales mentioned above.Good agreement with the characteristic time scales is indicated by observed times fitting to the dashed x = y -lines (Figure 10).
 In regime 1 (high Pe),  90% is proportional to   (Figure 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 fluidsolid heat transfer builds up and the solid temperature is still low.This corresponds to the time for traveling the full distance H.During stage 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   ∝   .For most cases it is controlled by diffusion through the solid over distances of order 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 time scale based on H (not shown).Inspecting this model shows that during stage 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 still non-equilibrium is present.
 In regime 2 (low Pe and  < 0.1, i.e.  > 10) the time for reaching ∆  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 time scale across the height of the system (Fig. 10b).
 In regime 3, (low Pe and high G (small H)), time for reaching ∆  is similar 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 order 1 and the advective timescale starts to take over.The time for reaching steady state (Fig. 10b, crosses) varies linearly with   ∝   .Clearly, it is also controlled by diffusion.

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 Peclet number , the initial thermal gradient G (= 1/H), and the melt fraction .Typical stages of melt flow for mid-ocean ridges include stage 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), stage b), merging melt channel or vein systems with high (> 10 -20%) porosity channels identified as dunite channels after complete melt extraction (Kelemen et al., 1997), and stage c), propagating dykes or other volcanic conduits.Let's assume typical overall melt fractions of 1% to 20% for stages b) and c).Schmeling et al. (2018)  As we preferably use the melt pore dimension   in our scalings (Eq.9a and 10a) we need to relate it to the solid phase dimension   by using Using ( 35), (9a), and ( 16) we arrive at the Peclet number used here m/s.With these paramterers   -numbers for the three stages can be estimated as 10 -7 to 10 -5 for stage a), 10 -5 to 10 -4 for stage b) at depths where channel distances are of order 0.1 m, and 10 -4 to 0.1 at shallower depths where the channel distances have increased to the order of 1m to 100 m, and >10 5 for the dyke stage c).
For mid ocean ridge settings we assume  of the order 1 to 10 km, and use Eq. ( 35) to insert typical   -values (increasing from 10 -4 m (stage a), interstitial melts) to 10 -3 m to 10 -1 m for the channeling stage b) (see Schmeling et al., 2018) to >10 m for the dyke stage c).The resulting Peclet number (Eq.37) is of the order 10 -3 to 0.5 for stage a), order 10 -2 during the early stage b) and order 10 -2 to 1 during the later stage b) appropriate for dunite systems, and order 10 4 to 10 7 for the dyke 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' With the derived estimates for the three stages, G' is of the order 10 -6 to 10 -2.7 for stage a), 10 -4 -10 -2.5 increasing to 10 -4 -0.6 for stage 2), and 10 -5 -10 -2 for the dyke stage c).These resulting stages for  and ′ are indicated in the regime diagram (Fig. 4).Starting from interstitial melts at full thermal equilibrium, channeling and veining may result in moderate thermal nonequilibrium at sufficiently high thermal gradients, while after transition to dyking 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, Peclet numbers are expected to be smaller, but non-dimensional thermal gradients (Eq.38) might be larger, resulting in a downward and rightward shift of the natural stages indicated in Figure 4.
To make our scaling laws and time scales for reaching maximum thermal non-equilibrium more accessible it is worth writing them in dimensional form.First, to estimate the Peclet number of a natural system combining Eq. ( 9) and ( 16) gives indicating that for very small or very large melt fractions Pe becomes very small.One may use Eq. ( 11) or ( 12) to write Pe also in terms of pore or grain dimensions df or ds, respectively.The scaling laws and characteristic time scales for the three regimes we found (Fig. 4) are in dimensional form:  Regime 1: For large Pe the maximum non-equilibrium temperature difference is simply equal to the imposed temperature difference, Δ  = Δ 0 , and the characteristic time to reach maximum non-equilibrium is simply  ℎ = / 0 , i.e. the total time of a fluid particle for passing through the system.
 Regime 2 and 3: For small Peclet number (  < √ √ 0 (1− 0 ) ) the maximum temperature difference scales like and the characteristic time for reaching this non-equilibrium scales with t0, i.e.
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.

Conclusions
In conclusion we showed that in magmatic systems characterized by two-phase flows of melts with respect to solid, thermalnon-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 Peclet number , the melt porosity , and the initial non-dimensional temperature gradient G in the system.The maximum possible non-equilibrium solidfluid temperature difference Δ  is controlled only by two nondimensional numbers: Pe and G.Both numerical and analytical solutions show that in a  −  -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   −   =  has been derived.
• In regime 3 (low Pe (<1) and G of order 1)) non-equilibrium scales with Pe and G and is depth-dependent, the scaling law is   −   =   (1 − ()) where M(z) depends on G.
Further conclusions include: • The time scales for reaching thermal non-equilibrium scale with the advective time-scale 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 Pe-G regime diagram reveals that a) interstitial melt flow is in thermal equilibrium, b) melt channeling as e.g.revealed by dunite channels may reach moderate thermal nonequilibrium, and c) the dyke regime is at full thermal non-equilibrium.
• In the studied setup G was constant leading to conservative estimates of thermal non-equilibrium.Any other depthdependent 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 watersandstone 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 nonequilibrium for anastomosing systems (Hart, 1993;Chevalier and Schmeling, in prep.).
) or (18) require two well defined boundary conditions each, Dirichlet (fixed temperature), Neumann (fixed thermal gradient), Robin (fixed sum of advective and conductive heat flux) 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.In an open outflow situation like our system neither the evolution of the temperature, the thermal gradient or the total (advective plus conductive) heat flux is known a priori, but depends 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.Only at later stages this boundary condition imposes artefacts in the temperatures field close to the top boundary.The limitations of this top boundary condition are tested and discussed in chapter 5.1.2.

Stage 1 :
and Model 3: H = 100, ϕ = 0.2.Figure2a and bshow 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), respectively.Figure2cshows the evolution of Tf and Ts with time at the top of the domain, for the same model 2 as in Figure2band for model 3 with a higher melt fraction ϕ = 0.2.Figure2dshows the evolution of (Tf -Ts) 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, evolve following three stages: During this transient stage the fluid temperature increases faster than the solid temperature (Fig.2a,b,c,e), and the temperature difference (Fig.2d,f) increases.During this stage, the fluid temperature increases rapidly at first, then the temperature increase slows down.As for the solid temperature, it first increases slowly, 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, thus the heat transfer term, increases too, discussed possible Peclet numbers for such systems based on a Darcy flow related Peclet number 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 100m -300 m for stage c) (dykes).Arguing for typical geometries, spreading rates and melt extraction ratesSchmeling et al. 2018 estimated the Darcy velocity lying between 10 -10 m/s and 10 -9

Figure 3 .
Figure 3. Maximum fluidsolid temperature differences   −   of numerical models (asterisks) with different parameters, plotted a) as a function of the Peclet number Pe for H = 10 and  = ., and b) as a function of the initial thermal gradient G for Pe = 1 and  = ..The solid lines give the analytic solutions.

Figure 4 .Figure 5 .
Figure 4. Main regimes of the maximum fluidsolid temperature differences   due to thermal non-equilibrium obtained by the analytical solution (Eq.22) in the parameter space of the Peclet number Pe and temperature gradient G.The asymptotic limits 645

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Temporal evolution of vertical profiles of (  −   ) for models with different Peclet numbers and model heights, i.e. different initial temperature gradients G = 1/H.In each panel the curves show (  −   )profiles for progressive times, the colors are cyclically varied with time from blue to yellow, starting with blue (bold curve).The first 5 curves of the Pe < 100 (respectively Pe = 100) models were taken with time increments of 1 (respectively 0.1), the later curves with 10 (respectively 1).The total time was 100 in all models with H = 10 and 500 in the models with H = 100.In each row the top boundary conditions is assumed as indicated 670

Figure 10 .
Figure 10.For evaluating time scales the numerically determined times of models with various parameters Pe and H 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 ∆  are plotted against either the advective time scale tadvH based on model height H for regime 1 models, or against the scaling time t0 for regime 2 models, or against the diffusive time scale tdiffH based on the model height H. b) times for reaching steady states are plotted against the characteristic diffusive time scales, tdiffH, based on model height H for all 3 regimes.Models close to the dashed line (y = x) are in best agreement with the characteristic times.In a) the Regime 2 times are taken dimensional by multiplying the observed times and the non-dimensional scaling time t0' = 1 by some arbitrary dimensional times t0.
To estimate Peclet numbers as defined here (Eq.36) typical interfacial thermal boundary layer thicknesses  are needed.As the thermal interfacial heat exchange intrinsically is time-dependent a good estimate is  =  ℎ √ 0  (in dimensional form) where  ℎ is a constant for a thermal boundary layer, equal to 2.32 for a cooling half space.Assuming that the characteristic time can be expressed by the (dimensional) fluid velocity  0 and system height , i.e. by  =   0 ⁄ = /  , we may express   in terms of the Peclet number PeD.With the resulting  and  we arrive at the following Peclet number (H and df are dimensional or non-dimensional):