A porous flow approach to model thermal non-equilibrium applicable to melt migration

A porous flow approach to model thermal non-equilibrium applicable to melt migration Abstract We develop an approach for heat exchange between a fluid and a solid phase of a porous medium where the temperatures of the fluid and matrix are not in thermal equilibrium. The formulation considers moving of the fluid within a resting or deforming porous matrix in an Eulerian coordinate system. The approach can be applied, for example, to partially molten systems or to brine transport in porous rocks. We start from an existing theory for heat exchange where the energy conservation equations for the fluid and the solid phases are separated and coupled by a heat exchange term. This term is extended to account for the full history of heat exchange. It depends on the microscopic geometry of the fluid phase. For the case of solid containing hot, fluid-filled channels, we derive an expression based on a time-dependent Fourier approach for periodic half-waves. On the macroscopic scale, the temporal evolution of the heat exchange leads to a convolution integral along the flow path of the solid, which simplifies considerably in case of a resting matrix. The evolution of the temperature in both phases with time is derived by inserting the heat exchange term into the energy equations. We explore the effects of thermal non-equilibrium between fluid and solid by considering simple cases with sudden temperature differences between fluid and solid as initial or boundary conditions, and by varying the fluid velocity with respect to the resting porous solid. Our results agree well with an analytical solution for non-moving fluid and solid. The temperature difference between solid and fluid depends on the Peclet number based on the Darcy velocity. For Peclet numbers larger than 1, the temperature difference after one diffusion time reaches 5 per cent of $$\tilde{T}$$ or more ($$\tilde{T}$$ is a scaling temperature, e.g. the initial temperature difference). Thus, our results imply that thermal non-equilibrium can play an important role for melt migration through partially molten systems where melt focuses into melt channels near the transition to melt ascent by dykes. Our method is based on solving the convolution integration for the heat exchange over the full flow history, which is numerically expensive. We tested to replace the heat exchange term by an instantaneous, approximate term. We found considerable errors on the short timescale, but a good agreement on the long timescale if appropriate parameters for the approximate terms are used. We derived these parameters which may be implemented in fully dynamical two-phase flow formulations of melt migration in the Earth. Hydrothermal systems, Mantle processes, Heat generation and transport, Mid-ocean ridge processes, Magma genesis and partial melting 1 INTRODUCTION Melt segregation from the deeper Earth to the upper crust or hydrothermal fluid flow in sedimentary basins are important examples of two-phase flows in geoscientific studies. Numerical modeling of such settings to quantify fluid flow and heat transport are commonly based on the Darcy equation describing fluid flow through a porous medium on a macroscopic scale (i.e. flow paths are not resolved). Porosity/permeability may be an intrinsic property of the rock or might be created by the formation of melt. Fluid flow is driven by thermal or compositional density differences between the fluid and the solid phases and a pressure gradient related to the state of stress of the fluid-filled rock. This typically applies to melt segregation within a partially molten rock volume, for example, beneath a mid-ocean ridge or within a plume. The general assumption for modeling fluid flow through a resting rock matrix (or through a matrix with a much slower velocity than the fluid) is thermal equilibrium, that is, temperature may vary spatially and temporally, but is equal in both phases. For hydrothermal circulation in the crust, the solid rock matrix is usually assumed to be undeformable and at rest. In contrast, for melt segregation in the Earth's crust or mantle the rock matrix is commonly described by a deformable viscous or viscoplastic material. In this case, one needs to address the fluid dynamics of a medium consisting of two immiscible fluids with different flow velocity fields. Such a problem is commonly referred to as two-phase flow. While the general theory for hydrothermal fluid flow in a resting matrix is well developed, the theory of two-phase flow is still in the focus of recent research. McKenzie (1984) and Spiegelman (1993) developed the theory of two-phase flow in the limit of high fluid-matrix viscosity contrasts between melt and matrix and low melt fractions. Bercovici et al. (2001) extended the theory to the case of large fluid fractions. Rudge et al. (2011) included chemical disequilibrium and melting. All these formulations assumed local thermal equilibrium based on the assumption that the thermal diffusion time on the grain scale, or more general, on the scale of a characteristic fluid phase length, is much shorter than the timescale of macroscopic flow and heat transport. However, there might exist cases in which the assumption of local thermal equilibrium is violated. Such cases might be rapid channeling of melt (e.g. Stevenson 1989; Golabek et al. 2008; Katz 2008), and the formation of melt-filled dykes in the crust (Rubin 1998). Also, water-filled fractures in the upper crust, important for hydrothermal fluid circulation, might be associated with diffusive length scales which violate the thermal equilibrium criteria. While thermal equilibrium on the scale of connected pores might be a valid assumption, it is certainty violated on the scale of fluid-filled dykes or fractures for non-isothermal systems. In such cases, hot fluids might intrude into the cooler rock matrix over large distances and exchange heat over some interval in time. Particularity, hot melt will intrude into a rock matrix even if the temperature of the rock matrix is below the melting temperature of the fluid material. The intrusion will continue as long as the temperature of the fluid is above its solidus. In the case of hydrothermal equilibrium fluid flow, the characteristic length of fluid transport in the medium can be related to the grain (or pore) size (Fig. 1a), in the case of non-hydrothermal equilibrium the characteristic length would be a characteristic host rock width ds (Fig. 1b), which is the averaged distance between adjacent solid–fluid interfaces through the solid phase. While the porosity φ might be the same in both cases, the conditions for heat transport differ. Figure 1. View largeDownload slide Sketch to illustrate the validity of the assumption of thermal equilibrium in porous flow: (a) fluid flow through pores and (b) fluid flow through dykes or fractures. Figure 1. View largeDownload slide Sketch to illustrate the validity of the assumption of thermal equilibrium in porous flow: (a) fluid flow through pores and (b) fluid flow through dykes or fractures. Amiri & Vafai (1994), Minkowycz et al. (1999), Nield & Bejan (2006), de Lemos (2016) and others presented approximate formulations for thermal non-equilibrium fluid flow considering only the actual, instantaneous temperature difference between fluid and solid phases for heat conduction between the phases and invoking appropriate heat transfer coefficients. In this paper, we follow a more fundamental physical approach and consider not only the actual temperature difference between the two phases, but include the time-dependent effects due to differential movement of the fluid and the solid (e.g. by melt migration) and related earlier temperature changes and heat exchange between the phases. In other words, we will explicitly account for the full time-dependent history of temperature differences between fluid and solid. This history is microscopically hidden within time-dependent thermal boundary layers with microscopically varying temperature gradients normal to the solid–fluid interfaces. Such gradients will be captured by our time-integration approach. Our theory is applicable to both, fluid flow and melt segregation, in a resting undeformable matrix, and to fluid flow or melt segregation in a deformable matrix. 2 OUTLINE OF THE THEORY In our general formulation of two-phase flow with local thermal non-equilibrium, we do not make a physical distinction between the two moving phases, thus the formulation will be phase symmetric. However, for semantic reasons we will call one phase ‘fluid’ with the subscript f, and the other phase ‘solid’ with the subscript s. The phases may differ in their physical properties, particularly in their densities, viscosities and thermal properties. The volume fraction occupied by the fluid will be addressed by the term ‘porosity’, though it is not necessarily a porosity in its strict sense. We distinguish between a microscopic and a macroscopic scale. On the microscopic scale, the two phases, fluid and solid, have distinct geometries and are separated by interfaces, and the temperature in each phase may vary spatially. However, in our approach, we do not treat this interface explicitly, but, on a macroscopic scale, we regard the two phases as continua which exchange heat (see Fig. 2). Thus, the temperature in each phase may vary spatially and may be different in each phase. On the macroscopic scale, the amount of heat exchange depends on the previous history of the temperature difference between solid and fluid as well as on variables such as the actual local temperature difference, thermal conductivities of the phases and a geometric factor characterizing the ratio of the interface surface to the volume. Figure 2. View largeDownload slide Sketch to illustrate the difference between (a) microscopic and (b) macroscopic scales: (a) distinct fluid-filled structures are resolved and the temperature in each phase might be different and (b) a larger region is considered, where fluid is only present in the central square box (grey pattern) with volumetric fraction φ. Both phases are considered as continua which locally can exchange heat and where the temperatures in the phases can be different. Figure 2. View largeDownload slide Sketch to illustrate the difference between (a) microscopic and (b) macroscopic scales: (a) distinct fluid-filled structures are resolved and the temperature in each phase might be different and (b) a larger region is considered, where fluid is only present in the central square box (grey pattern) with volumetric fraction φ. Both phases are considered as continua which locally can exchange heat and where the temperatures in the phases can be different. More specifically, our approach makes use of a time-dependent reference heat exchange function between the solid and fluid ($$q_s^*( t )$$, see below), which has to be determined on the microscopic scale prior to solving the full problem. When solving the full problem on the macroscopic scale by separating the heat equation into one for the solid and one for the fluid, the heat exchange between the two phases will take into account the actual temperature difference between the phases, as well as all previous temperature differences. These have to be convoluted with the reference heat exchange function to account for thermal relaxation. Thus, our approach requires storing the actual and all previous temperatures of both phases for the complete run. Our approach contains some simplifications compared to solving the full microscopic advection-diffusion problem for porous flow: the detailed 3-D microscopic solid and fluid geometries are replaced by equivalent geometries such as layers, spheres or infinite half-spaces. Thus, the detailed geometry of a porous network cannot be captured. No local heat advection in the fluid by differential flow on the pore scale, only the net Darcy flow is accounted for. We assume that the two phases obey the equations of mass and energy conservation (see e.g. Schmeling 2000). In the present approach, we neglect the melting and solidification process, however, this process might be introduced later by invoking the Stefan problem. The velocities of the solid and the fluid phases can be the results of solving the momentum equations for these phases, for example, a Stokes equation for the deformable porous solid and a Darcy equation for the fluid, but the way they are obtained is not of concern for this study. 2.1 Basic equations To understand how heat is exchanged between the fluid and the solid phases, we will first examine the process on the microscopic scale (see Fig. 3a). We assume that at a certain time the locally averaged temperature of the fluid phase in a medium of porosity φ (i.e. the volumetric fraction of the fluid) is given by Tf, that of the solid is Ts, and that of the interface between fluid and solid is Tb. Local averaging is meant as averaging the temperatures over the respective fluid and solid volumes, and the interface, on a length scale larger than the characteristic pore scale (ds in Figs 1 and 2), but smaller than the macroscopic length scale (L in Fig. 2). Conservation of energy of the matrix phase is given by (e.g. Nield & Bejan 2006; de Lemos 2016)   \begin{equation}\left( {1 - \varphi } \right){\rho _s}{c_{ps}}\left( {\frac{{\partial {T_s}}}{{\partial t}} + {{\vec{v}}_s} \cdot \vec{\nabla }{T_s}} \right) = \vec{\nabla }{\rm{\ }} \cdot \left( {\left( {1 - \varphi } \right){\lambda _s}\vec{\nabla }{T_s}} \right) + {Q_b},\end{equation} (1) and that of the fluid phase is   \begin{equation}\varphi {\rho _f}{c_{pf}} \left( {\frac{{\partial {T_f}}}{{\partial t}} + {{\vec{v}}_f} \cdot \vec{\nabla }{T_f}} \right) = \vec{\nabla } \cdot \left( {\varphi {\lambda _f}\vec{\nabla }{T_f}} \right) - {Q_b},\end{equation} (2) Figure 3. View largeDownload slide (a) Microscopic sketch of the thermal non-equilibrium two-phase flow. The temperatures Ts, Tf, Tband the velocities are variable on this microscopic scale, but constant on the macroscopic scale. (b) Temporal evolution of the local temperature field (shown in the grey ellipsis in Fig. 3a)as due to a change in the macroscopic temperature difference between solid and fluid, as a function of the distance x from the interface. Thick black line for time t = 0 and thinner black lines for successive times. Figure 3. View largeDownload slide (a) Microscopic sketch of the thermal non-equilibrium two-phase flow. The temperatures Ts, Tf, Tband the velocities are variable on this microscopic scale, but constant on the macroscopic scale. (b) Temporal evolution of the local temperature field (shown in the grey ellipsis in Fig. 3a)as due to a change in the macroscopic temperature difference between solid and fluid, as a function of the distance x from the interface. Thick black line for time t = 0 and thinner black lines for successive times. where Qb is the interphase heat exchange (IPHE) rate (i.e. the thermal coupling term of the fluid and solid phase across the interfaces) with the dimension Jm−3 s−1; all other variables are explained in Table 1. In our formulation of eqs (1) and (2), we neglect heat generation due to internal heating, latent heat and dissipation (but see Rudge et al. (2011) for a more complete formulation for the thermal equilibrium case). Note, that in non-equilibrium two-phase flow, we need two separate equations for heat transport. Table 1. Symbols, their definitions and physical units used in this study. Symbol  Definition  Units  as  Geometry factor depending on the dimension of the problem  –  cps  Specific heat capacity of the solid phase  J kg−1 K−1  cpf  Specific heat capacity of the fluid phase  J kg−1 K−1  d  Grain size  M  df  Characteristic length of the fluid phase  M  ds  Characteristic length of the solid phase  M  L  Macroscopic length scale  M  l0  Wavelength for fracture/dyke spacing  M  Pe  Peclet number  –  PeD  Darcy flow related Peclet number  –  PeL  Macroscopic Peclet number  –  Qb  Interphase heat exchange rate (IPHE)  J s−1 m−3  Δq  Heat flow through the fluid–solid interface  W m−2  qs  Specific heat flow from the fluid to the solid phase, normal to the interface  W m−2  $$q_s^*$$  Specific heat flow response function to a temperature step function  W m−2 K−1  qf  Specific heat flow from the solid to the fluid phase, normal to the interface  W m−2  $${q_{{s_{{\rm{inst}}}}}}$$  Instantaneous specific heat flow from the solid to the fluid phase  W m−2  S  Geometry term, describing the interface area per volume of the solid  m−1  Sb  Surface area of all fluid–solid interfaces in a volume  m2  T  Temperature  °C  Tb, Tb0  Temperature and initial temperature of the fluid–solid interface  °C  Tf, Tf0  Temperature and initial temperature of the fluid phase  °C  Ts, Ts0  Temperature and initial temperature of the solid phase  °C  t, t0  Time and initial time  S  tov  ‘Overtake time’ when advective propagation prevails conductive  S  Vg  Volume containing fluid-filled fractures/dykes (entire volume)  m3  vD  Darcy velocity  m s−1  Vs  Volume of the solid phase  m3  $${\vec{v}_f}$$  Velocity of the fluid phase  m s−1  $${\vec{v}_s}$$  Velocity of the solid phase  m s−1  vdiff  Velocity of conductively moving temperature front  m s−1  W  Heat content per volume  J m−3  Ws, Wf  Heat content per volume of the solid phase and of the fluid phase  J m−3  $$\vec{x}$$  Position vector  M  x  Space coordinate normal to the fractures/dykes  M  xinter  x position of the fluid–solid interface  M  z  Space coordinate in direction of the fractures/dykes  M  δM  thickness of the thermal boundary layer  M  δT  Small temperature increment  K  ΔT0  Initial temperature difference between the fluid and the solid phase  K  κs, κf  Thermal diffusivity of the solid phase and of the fluid phase  m2 s−1  λs, λf  Thermal conductivity of the solid phase and the fluid phase  W m−1 K−1  $$\bar{\lambda }$$  Mean thermal conductivity  W m−1 K−1  ρs, ρf  Density of the solid phase and of the fluid phase  kg m−3  φ  Porosity  –  $${\vec{\psi }_s}$$  Particle path in the solid phase  M  $${\vec{\psi }_{s0}}$$  Initial starting point at the flow path  M  $${\vec{\psi }_f}$$  Particle path in the fluid phase  M  Symbol  Definition  Units  as  Geometry factor depending on the dimension of the problem  –  cps  Specific heat capacity of the solid phase  J kg−1 K−1  cpf  Specific heat capacity of the fluid phase  J kg−1 K−1  d  Grain size  M  df  Characteristic length of the fluid phase  M  ds  Characteristic length of the solid phase  M  L  Macroscopic length scale  M  l0  Wavelength for fracture/dyke spacing  M  Pe  Peclet number  –  PeD  Darcy flow related Peclet number  –  PeL  Macroscopic Peclet number  –  Qb  Interphase heat exchange rate (IPHE)  J s−1 m−3  Δq  Heat flow through the fluid–solid interface  W m−2  qs  Specific heat flow from the fluid to the solid phase, normal to the interface  W m−2  $$q_s^*$$  Specific heat flow response function to a temperature step function  W m−2 K−1  qf  Specific heat flow from the solid to the fluid phase, normal to the interface  W m−2  $${q_{{s_{{\rm{inst}}}}}}$$  Instantaneous specific heat flow from the solid to the fluid phase  W m−2  S  Geometry term, describing the interface area per volume of the solid  m−1  Sb  Surface area of all fluid–solid interfaces in a volume  m2  T  Temperature  °C  Tb, Tb0  Temperature and initial temperature of the fluid–solid interface  °C  Tf, Tf0  Temperature and initial temperature of the fluid phase  °C  Ts, Ts0  Temperature and initial temperature of the solid phase  °C  t, t0  Time and initial time  S  tov  ‘Overtake time’ when advective propagation prevails conductive  S  Vg  Volume containing fluid-filled fractures/dykes (entire volume)  m3  vD  Darcy velocity  m s−1  Vs  Volume of the solid phase  m3  $${\vec{v}_f}$$  Velocity of the fluid phase  m s−1  $${\vec{v}_s}$$  Velocity of the solid phase  m s−1  vdiff  Velocity of conductively moving temperature front  m s−1  W  Heat content per volume  J m−3  Ws, Wf  Heat content per volume of the solid phase and of the fluid phase  J m−3  $$\vec{x}$$  Position vector  M  x  Space coordinate normal to the fractures/dykes  M  xinter  x position of the fluid–solid interface  M  z  Space coordinate in direction of the fractures/dykes  M  δM  thickness of the thermal boundary layer  M  δT  Small temperature increment  K  ΔT0  Initial temperature difference between the fluid and the solid phase  K  κs, κf  Thermal diffusivity of the solid phase and of the fluid phase  m2 s−1  λs, λf  Thermal conductivity of the solid phase and the fluid phase  W m−1 K−1  $$\bar{\lambda }$$  Mean thermal conductivity  W m−1 K−1  ρs, ρf  Density of the solid phase and of the fluid phase  kg m−3  φ  Porosity  –  $${\vec{\psi }_s}$$  Particle path in the solid phase  M  $${\vec{\psi }_{s0}}$$  Initial starting point at the flow path  M  $${\vec{\psi }_f}$$  Particle path in the fluid phase  M  View Large The main problem is the derivation of the IPHE term. To determine an expression for the IPHE, we start on the microscopic scale by considering the time-dependent evolution of the temperature field on either side of a fluid–solid interface with the (hot) fluid phase moving with a differential velocity compared to the solid phase (see Fig. 3a). The temperature difference between solid and fluid may also be caused by different conductive heat fluxes through the two phases due to different thermal parameters. The principal assumption is that at a certain time t0 an initial temperature difference ΔT exists between the solid and fluid phases (e.g. due to hot melt intrusion). This temperature difference initiates a temporally and spatially varying temperature field in either phase on the microscopic scale. For successive times (t > t0), the initial temperature difference will equilibrate (see Fig. 3b). Thus, any temperature difference between a solid and a fluid phase generally represents a transient stage. If we assume a positive temperature change in a hot fluid or melt, it acts as a heat source in the solid (eq. 1) and a sink in the fluid (eq. 2, i.e. the hot fluid is losing heat to the colder solid). It depends on the local heat flux across the interface and the specific interfacial area S of dimension m−1 describing the interface area per volume   \begin{equation}{Q_b} \left( {\vec{x},t} \right) = \ S \cdot {q_s}\left( {\vec{x},t} \right).\end{equation} (3) qs is the heat flow normal to the interface from the fluid phase into the solid phase (see Fig. 3a). Since we only consider the normal component of the heat flow it is treated in the following as a scalar variable. Qb could be determined in a similar way by using the interface normal of the heat flow from the solid phase to the fluid phase qf. We will first derive an expression for $${q_s}( {\vec{x},t} )$$ and further down discuss the specific interfacial area term S. 2.2 Interface heat flux In the following, we treat the heat flux resulting from a sudden temperature change as described by the ‘heat flow response function to a temperature step function’, denoted by an asterisk $$q_s^*( t )$$  \begin{equation}{q_s} \left( t \right) = \ \Delta {T_0} \cdot q_s^*\left( t \right).\end{equation} (4) The expression for $$q_s^*( t )$$ depends on the application, i.e. the choice of the phase geometries. A possible choice could be the assumption of a cooling half-space ($$q_s^*( t ) = {\lambda _s}\frac{1}{{\sqrt {\pi {\kappa _s}t} }}$$, e.g. Turcotte & Schubert 2002). Another possibility is the assumption of a series of cooling (vertical or horizontal) sheets (for which we will derive the expression further down). As a result of an initial temperature difference  ΔT0 at all interfaces in a volume, a subsequent heat flux from the fluid into the solid,$${\rm{\ \Delta }}{T_0}q_{s0}^*( {\tilde{t}} )$$, is generated, where $$\tilde{t}$$ is the subsequent time starting with t0, that is, $$\tilde{t} = \ t - {t_0}$$. For subsequent times  t > t0, we determine the heat flux into the solid by considering any temporal changes δT(τ) of the macroscopic fluid–solid temperature difference along a solid volume flow path (see Fig. 4) taking place at any time in the history of the solid particle τ < t. Each of these incremental temperature changes cause an additional time-dependent heat flow into the solid, which, at the time t, has decayed to $$d{q_{s,\tau }}( t ) = \ \delta T( \tau )q_s^*( {t - \tau } )$$. The full interface heat flux at time t can thus be determined by summing up $${\rm{\Delta }}{T_0}q_{s0}^*( {\tilde{t}} )$$ and an integral over all dqs, τ(t) having emerged at all times τ < t. The problem is that this integration has to take into account the history of solid and fluid flow paths, and that the temperature changes δT(τ) need some further attention. Figure 4. View largeDownload slide Solid volume flow path $${\vec{\psi }_s}( \tau )$$and fluid volume flow paths $${\vec{\psi }_{fi}}( \tau )$$, where the fluid flow paths cross the solid flow path at different times. The fluid flow path $${\vec{\psi }_{fn}}( \tau )$$merges the solid path at $$( {\vec{x},t} )$$. At previous time τi , the solid particle has been exposed to different fluid particles, denoted by the fluid flow paths. At those times, the solid particle has experienced incremental temperature changes $${ {\delta T} |_{{{\vec{\psi }}_s}( \tau )}}$$ whose subsequent response needs to be determined along the solid flow path all the way to $$( {\vec{x},t} )$$. Figure 4. View largeDownload slide Solid volume flow path $${\vec{\psi }_s}( \tau )$$and fluid volume flow paths $${\vec{\psi }_{fi}}( \tau )$$, where the fluid flow paths cross the solid flow path at different times. The fluid flow path $${\vec{\psi }_{fn}}( \tau )$$merges the solid path at $$( {\vec{x},t} )$$. At previous time τi , the solid particle has been exposed to different fluid particles, denoted by the fluid flow paths. At those times, the solid particle has experienced incremental temperature changes $${ {\delta T} |_{{{\vec{\psi }}_s}( \tau )}}$$ whose subsequent response needs to be determined along the solid flow path all the way to $$( {\vec{x},t} )$$. We consider that the fluid and possibly also the solid move on independent paths (in a fixed Eulerian coordinate system). The heat fluxes and temperatures discussed above, however, are related to the heat contents of fluid and solid volumes, which will be carried away by either flow. The position of a moving solid volume particle can be described by a flow path (time-dependent position vector) $${\vec{\psi }_s}( t )$$ (see Fig. 4), starting at an initial position $${\vec{\psi }_{s0}}$$. We have to consider all incremental fluid–solid temperature changes (δT = δTf − δTs) as incremental step functions along such a solid particle volume flow path, which are not associated with the instantaneous interphase heat exchange (IPHE), i.e. the heat exchange along the solid flow line at a particular time (denoted further down by τ). In other words, δTs and δTf may be due to some ‘external’ macroscopic heat transport effects such as advection, conduction, heat generation and phase changes. In particular, we need to consider processes which act differently in the two phases, for example, advection with different $${\vec{v}_s}$$ and $${\vec{v}_f}$$, conduction with different thermal conductivities. For the incremental temperature change δT along the path at times τ ≤ t, we can write   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} - \left( {{{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}}} \right)} \right)\ \delta \tau ,\end{equation} (5) where $${ {\frac{{d{T_f}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ are the entire temperature changes in the solid and the fluid phases along the path $${\vec{\psi }_s}( \tau )$$ (i.e. following the solid particle), respectively, and $${ {\frac{{d{T_f}}}{{dt}}} |_{{Q_b},{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{Q_b},{{\vec{\psi }}_s}( \tau )}}$$ are the changes due to the instantaneous interphase heat exchange. While $${ {\frac{{d{T_f}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ can simply be evaluated along the flow path, an expression for $$( { - {{ {\frac{{d{T_f}}}{{dt}}} |}_{{Q_b},{{\vec{\psi }}_s}( \tau )}} + {{ {\frac{{d{T_s}}}{{dt}}} |}_{{Q_b},{{\vec{\psi }}_s}( \tau )}}} )$$ needs to be derived. The heat contents W per volume (including both phases) at any instant of time is given for the solid and the fluid phases by   \begin{equation}{W_s} = \left( {1 - \varphi } \right)\ {\rho _s}{c_{Ps}}{T_s} \,{\rm{and}}\,{W_f} = \varphi {\rho _f}\ {c_{Pf}}{T_f}.\end{equation} (6) These heat contents change locally with time, and the contribution due to the instantaneous interphase heat exchange rate (denoted by subscript Qb) can written in terms of the IPHE (eq. 3) as   \begin{equation}S \cdot \ {q_s} = {\left. {\frac{{d{W_s}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = \ - \ {\left. {\frac{{d{W_f}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {1 - \varphi } \right)\ {\rho _s}{c_{Ps}}\ {\left. {\frac{{d{T_s}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = - \varphi {\rho _f}{c_{Pf}}{\left. {\frac{{d{T_f}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}}.\end{equation} (7) Rearranging for the temperature change rates and inserting them in eq. (5) leads to   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {\frac{{S \cdot {q_s}}}{{\varphi {\rho _f}{c_{Pf}}}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {\frac{{S \cdot {q_s}}}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}}} \right)\delta \tau \,,\end{equation} (8) using again eq. (3), this is   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {{Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}}} \right)\ \delta \tau .\end{equation} (9) To consider all incremental temperature changes δT along a solid volume path, we sum up the initial heat flow response, $${\rm{\Delta }}{T_0}q_{s0}^*( {t - {t_0}} )$$, and the $$d{q_{s,\tau }}( t ) = \ \delta T( \tau )q_s^*( {t - \tau } )$$ increments integrated over all previous times τ. This integral, $$\mathop \smallint \nolimits_{{t_0}}^t q_s^*( {t - \tau } ){ {\frac{{\delta T}}{{\delta \tau }}} |_{{{\vec{\psi }}_s}( \tau )}}\,d\tau $$, has to account for the temperature changes δT(τ) along the solid volume path $${\vec{\psi }_s}( \tau )$$ from the initial time t0 to the actual time t. Taking   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = {\left. {\frac{{\partial T}}{{\partial t}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot {\left. {\frac{{dt}}{{d{{\vec{\psi }}_s}}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot d{\psi _s} = {\left. {\frac{{\partial T}}{{\partial t}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot {\left. {\frac{1}{{{{\vec{v}}_s}}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot d{\psi _s},\end{equation} (10) where $${\vec{v}_s}$$ is the velocity of the solid phase, and using eq. (9), the sum of the initial heat flow response and the integrated heat flow increments yields   \begin{equation}{q_s}\!\left( {{{\vec{\psi }}_s}\left( t \right)} \right) = \ \Delta {T_0}\ q_{s0}^*{\left( {t - {t_0}} \right)_{{{\vec{\psi }}_{s0}}}} + \int_{{{{\vec{\psi }}_{s0}}}}^{{{{\vec{\psi }}_s}\left( t \right)}}{{q_s^*\left( {t - \tau \left( {\psi _s^\prime} \right)} \right)}} \cdot \ {\left( {\frac{{d{T_f}}}{{dt}} - \frac{{d{T_s}}}{{dt}} + {Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right)_{\psi _s^{'}}} \cdot {\left. {\frac{1}{{{{\vec{v}}_s}}}} \right|_{\psi _s^{'}}}d\psi _s^{'}.\end{equation} (11)$$\tau ( {{{\vec{\psi }}_s}} )$$ is the time at which the considered solid particle volume passed a previous position $${\vec{\psi }_s}$$. In case of a resting solid and only a moving fluid, eq. (11) degenerates to a convolution integral   \begin{equation}{q_s}\!\left( {\vec{x},t} \right) = \ \Delta {T_0}\ q_{s0}^*\left( {t - {t_0}} \right) + \int_{{{t_0}}}^{t}{{q_s^*\left( {t - \tau } \right)}} \cdot \ {\left( {\frac{{d{T_f}}}{{dt}} - \frac{{d{T_s}}}{{dt}} + {Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right)_\tau }d\tau ,\end{equation} (12) eq. (12) allows determining the interphase heat exchange rate at a given position and time by integrating over all previous temperature changes of the fluid, which has passed that position, as well as integrating the history of solid temperature changes and heat exchange at that position. 2.3 Heat exchange for a medium containing fluid layers When discussing eq. (4), we already mentioned the choice of a cooling half-space solution for $$q_s^*( t )$$. This is a suitable choice if the underlying assumption is related to widely spaced fluid channels, dykes or fractures (or fluid layers in general) with heat exchange to the solid but no heat exchange in between the fractures. For densely spaced fluid layers such an assumption does not hold. Here, we develop an expression for heat flow from the fluid into the solid for a case of densely spaced fluid layers (or vertical dykes) where the layer spacing influences the overall heat exchange. If we recall Fig. 3(b), a case of a swarm of large vertical fluid layers can be constructed as follows: since the layers are large in vertical direction and we regard only the normal component of the heat flow from the fluid-filled layers into the solid, we can treat the problem as 1D. For a swarm of fluid layers of equal spacing ds and width df, the initial or any incremental temperature difference between fluid and solid which is advected to a certain macroscopic position can be regarded as an initial box-shaped temperature function on microscopical scale. The related initial temperature as a function of a space coordinate x normal to the fluid layer is described by the Fourier series of a rectangular half-wave with an initial temperature difference between fluid and solid phases of ΔT0 and fluid layer spacing given by a wavelength l0 (see Fig. 5)   \begin{equation}{T_0} \left( x \right) = \frac{{{d_f} \cdot \Delta {T_0}}}{{{l_0}}}\ + \mathop \sum \limits_{n\ = \ 1}^\infty {\frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \frac{{{d_f}n\pi }}{{{l_0}}}\cos \frac{{2\pi nx}}{{{l_0}}}} .\end{equation} (13) Figure 5. View largeDownload slide Sketch to illustrate the assumption of a periodic temperature perturbation to model narrowly spaced fluid layers such as dykes or fractures filled with a hot fluid. Figure 5. View largeDownload slide Sketch to illustrate the assumption of a periodic temperature perturbation to model narrowly spaced fluid layers such as dykes or fractures filled with a hot fluid. Since the macroscopic porosity can be expressed by $$\varphi \ = \frac{{{d_f}}}{{{l_0}}}$$, eq. (13) is rewritten to   \begin{equation}{T_0} \left( x \right) = \ \Delta {T_0}\varphi + \mathop \sum \limits_{n\ = \ 1}^\infty {\frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \left( {n\pi \varphi } \right)\cos \frac{{2\pi n\varphi x}}{{{d_f}}}.} \end{equation} (14) Assuming the same thermal properties in the solid and fluid, the temperature for successive times t is determined by solving the 1D conductive heat equation $$\frac{{\partial T}}{{\partial t}} = {\kappa _s}\frac{{{\partial ^2}T}}{{\partial {x^2}}}$$ with the initial condition (eq. 14) for each Fourier term separately. The result is Tn (x, t) = T0n (x) · exp ( − bnt) where T0n(x) is the Fourier terms of the initial condition in eq. (14) and the coefficients bn are obtained as $${b_n} = {\kappa _s}{( {\frac{{2\pi n\varphi }}{{{d_f}}}} )^2}$$. Thus, the temperature in a continuum containing initially rectangular half-wave temperature differences is   \begin{equation}T ( {x,t} ) = \Delta {T_0}\varphi + \mathop \sum \limits_{n\ = \ 1}^\infty \frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \left( {n\pi \varphi } \right)\cos \frac{{2\pi n\varphi x}}{{{d_f}}}\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (15) Applying Fourier's law of heat conduction in one dimension $$q\ = \ - {\lambda _s}\frac{{dT}}{{dx}}$$ to eq. (15) gives the heat flow in such a continuum   \begin{equation}q ( {x,t} ) = {\lambda _s}\ \frac{{4\Delta {T_0}\varphi }}{{{d_f}}}\mathop \sum \limits_{n\ = \ 1}^\infty \sin \left( {n\pi \varphi } \right)\sin \frac{{2\pi n\varphi x}}{{{d_f}}}\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (16) Since we are interested in the heat flow qs from the fluid phase into the solid phase at the interface of the two phases, which is defined by $$x\ = {x_{{\rm{inter}}}}\ = \frac{{{d_f}}}{2}$$ (see Fig. 5), we evaluate eq. (16) at x = xinter  \begin{equation}{q_{s\_{\rm{Layer}}}}\ \left( t \right) = {\lambda _s}\ \frac{{4\Delta {T_0}\varphi }}{{{d_f}}}\mathop \sum \limits_{n\ = \ 1}^\infty {\sin ^2}\left( {n\pi \varphi } \right)\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (17) This is the exchange heat flux of the layered solid-fluid case. Eq. (17) can be written in the form of   \begin{equation}q_s^* \left(t \right) = \frac{1}{{\Delta {T_0}}}\ {q_{s\_{\rm{Layer}}}}\left( t \right)\end{equation} (18) (see eq. 4) to be used as the ‘heat flow response function to a temperature step function’ in eq. (11) and in eq. (12). 2.4 The geometry term The geometry term S in eq. (3) has the dimension m−1 and is defined by the surface area Sb of all fluid–solid interfaces in a volume element to the entire volume Vg  \begin{equation}S\ = \frac{{{S_b}}}{{{V_g}}}.\end{equation} (19) For spherical solid inclusions of diameter ϑ (i.e. radius ϑ/2), the interface surface is given by Sb = πϑ2 and the solid phase volume by $${V_s} = \frac{1}{6}\ \pi {\vartheta ^3}$$. Since Vg = Vs(1 − φ)− 1, S for spherical solid inclusions is S = 6ϑ− 1(1 − φ). One can do this little exercise for other geometries to find the general relation   \begin{equation}S\ = {a_s}\ d_s^{ - 1}\left( {1 - \varphi } \right),\end{equation} (20) where, more generally, ds is a scaling length of the solid phase of the composite material (i.e. the average distance between adjacent solid–fluid interfaces through the solid phase) and as ≈ 6 for 3-D (e.g. spherical solid volumes), as ≈ 4 for 2-D (e.g. cylindrical solid volumes) and as ≈ 2 for 1-D (e.g. solid and fluid layers). Our simplifying assumption for deriving eq. (20) implies complete wetness, that is, solid grains (or inclusions) are totally surrounded by fluid films. For the 2-D case of parallel fluid dykes, which we study here, this assumption is approximately valid. For other melt geometries such as melt pockets or tubes at the edges of solid grains, the wetness is less than 1 and a different formula than eq. (20) would emerge. 2.5 Non-dimensionalization of the equations For using the equations above in a numerical code and studying the principal behaviour of non-equilibrium two-phase flow, it is convenient to use the equations in non-dimensionalized form. We do this by choosing a set of (dimensional) scaling parameters ($$\skew4\tilde{d},\ \skew4\tilde{\kappa },\ \skew4\tilde{\rho },\ {\skew4\tilde{c}_p},\ \skew4\tilde{T}$$) and setting   \begin{equation}x = x{'} \cdot \tilde{d},\qquad t = t{'} \cdot \frac{{{{\tilde{d}}^2}}}{{\tilde{\kappa }}},\qquad \vec{v} = \vec{v}{'} \cdot \frac{{\tilde{\kappa }}}{{\tilde{d}}},\qquad T = T{'} \cdot \tilde{T},\qquad \lambda = \lambda {'} \cdot \tilde{\kappa }\,\tilde{\rho }\,{\tilde{c}_p},\end{equation} (21) where the apostrophe indicates non-dimensional variables. By specifying $$\tilde{d} = {d_s}\ ,\ \ \widetilde {\kappa \ } = {\kappa _s}\ ,\ \ \tilde{\rho } = {\rho _s}\ ,\ \ {\tilde{c}_p} = {c_{ps}}\ ,\ \ \tilde{T} = \ {\rm{\Delta \ }}{T_0} = {T_{f0}}\ - {T_{s0}}$$, we scale the initial temperature difference for the fluid phase and the solid phase to 1 and obtain for the heat flow, the specific interface S, and Qb  \begin{equation}q = {\kappa _s}\,{\rho _s}\,{c_{ps}}\frac{{\Delta {T_0}}}{{{d_s}}}q{'},\quad S = \frac{1}{{{d_s}}}S{'},\quad {Q_b} = {\kappa _s}\,{\rho _s}\,{c_{ps}}\frac{{\Delta {T_0}}}{{{d_s}^2}}{Q{'}_b}.\end{equation} (22) Non-dimensionalizing of eq. (18) and considering that $${d_f} = \frac{\varphi }{{( {1 - \varphi } )}}{d_s}$$ leads to   \begin{equation}q{_s^{*\prime} } \left( {t{'}} \right) = 4\left( {1 - \varphi } \right)\mathop \sum \limits_{n\ = \ 1}^\infty {\sin ^2}\left( {n\pi \varphi } \right)\exp \left( { - t{'}{{\left( {2\pi n\left( {1 - \varphi } \right)} \right)}^2}} \right).\end{equation} (23) We also need non-dimensionalized versions of eq. (1) for the solid phase   \begin{equation}\left( {1 - \varphi } \right)\ \left( {\frac{{\partial {T_s}^\prime }}{{\partial t{'}}} + {{\vec{v}}_s}{'} \cdot \vec{\nabla }{'}{T_s}{'}} \right) = \vec{\nabla }{'} \cdot \left( {\left( {1 - \varphi } \right)\vec{\nabla }{'}{T_s}{'}} \right) + {Q_b}^\prime ,\end{equation} (24) and eq. (2) for the fluid phase   \begin{equation}\varphi \ \left( {\frac{{\partial {T_f}{'}}}{{\partial t{'}}} + \vec{v}_f^{'} \cdot \vec{\nabla }{'}{T_f}{'}} \right) = \vec{\nabla }{'} \cdot \left( {\varphi \vec{\nabla }{'}{T_f}^\prime } \right) - {Q_b}^\prime .\end{equation} (25) The non-dimensionalized form of the geometry term is S' = as (1 − φ), and eq. (12) with the convolution integral is given by   \begin{equation}{q_s}^\prime \ \left( {\vec{x}{'},t{'}} \right) = q{_{s0}^{*\prime} }\left( {t{'} - {{t{'}}_0}} \right) + \mathop \smallint \nolimits_{{t_0}{'}}^{t{'}} q_s^*{'}\left( {t{'} - \tau {'}} \right) \cdot \ {\left( {\frac{{d{T_f}{'}}}{{dt{'}}} - \frac{{d{T_s}{'}}}{{dt{'}}} + {Q_b}{'}\left( {\frac{1}{{\varphi {\rho _f}{'}{c_{Pf}}{'}}} + \frac{1}{{\left( {1 - \varphi } \right)}}} \right)} \right)_{\tau {'}}}d\tau {'}.\end{equation} (26) 3. APPLICATION OF THE THEORY First, we determined the maximum Fourier coefficient nmax needed in the sum in eq. (23) to get a sufficient approximation for the response function for the case of parallel fluid layers (cf. Fig. 5). We tested values between 10 and 500 and found that >100 is necessary for reliable results. In this study, we used nmax = 500. Fig. 6 shows the non-dimensional heat flow response function $$q_s^{*{'}}$$ as a function of time and porosity in a log–log plot. Typically, three or two stages can be identified for all curves in Fig. 6, characterized by different slopes in the $$\log q_s^{*{'}} - \log t{'}$$ plot. The first stage, independent of the porosity, has a slope of $${\sqrt {t{'}} ^{ - 1}}$$ and represents the solution of two half-spaces of different initial temperature brought into sudden contact. Heat exchange is not yet influenced by the neighbouring fluid layer and the central regions in the solid and fluid layers are still at their initial temperatures. The second stage has a slope of about (t')− 3/2 (this can be verified by choosing porosities closer to 0 or 1) and reflects the stage in which the thermal boundary layers in the thinner layer (fluid or solid, depending on φ) have merged in the centre of the thinner layer, while the centre of the thicker layer is still not affected by the cooling/heating. In contrast to the first stage, the temperature at the boundary between the fluid and solid drops (or rises) with time, which leads to a stronger decay than proportional to $${\sqrt {t{'}} ^{ - 1}}$$. The last stage has a slope according to e− t'and describes the situation when heat flux rapidly ceases. For equally thick solid and fluid layers (φ = 0.5), the second stage is skipped, because the thermal boundary layers merge in the centres of both layers simultaneously. The differences of the curves in respect to porosity is caused by the difference in total heat contents of the fluid layers, since their volume is a function of the porosity while the heat exchange boundary is not in this 1-D case. Figure 6. View largeDownload slide The non-dimensional fluid–solid interface ‘heat flow response function to a temperature step function’ $$q_s^{*{'}}$$ (eq. 23) as a function of the non-dimensional time t'and the volumetric fraction of the fluid φ (‘porosity’). The sum in eq. (23) was truncated at used nmax = 500. The three different stages which are discussed in the text are indicated. Figure 6. View largeDownload slide The non-dimensional fluid–solid interface ‘heat flow response function to a temperature step function’ $$q_s^{*{'}}$$ (eq. 23) as a function of the non-dimensional time t'and the volumetric fraction of the fluid φ (‘porosity’). The sum in eq. (23) was truncated at used nmax = 500. The three different stages which are discussed in the text are indicated. To apply our theory for non-equilibrium two-phase flow, we studied two simple 1-D test cases. For all cases, we used $$q_s^{*{'}}( {t{'}} )$$ as given in eq. (23) and shown in Fig. 6, and solved eq. (26) for the heat flow from the fluid phase to the solid phase. With the result for $$q_s^{'}( {\vec{x},t} )$$, we determined the IPHE term $$Q_b^{'}( {\vec{x},t} )$$ (using S' and $$q_s^{'}( {\vec{x},t} )$$ in eq. 3) and setting as = 2. With this $$Q_b^{'}( {\vec{x},t} )$$, we solved the coupled heat transport equations (eqs 24 and 25) for the fluid and the solid phases, respectively. We studied the following test cases and assumed for the sake of simplicity the same thermal conductivity for both phases, λf = λs: Case 1: solid and fluid phases in rest. This is accomplished by setting $${\vec{v}_f} = \ 0$$ in eq. (25). Case 2: solid phase in rest, fluid phase moving. Case 2a evaluates the full convolution integral $$\mathop \smallint \nolimits_{{t_0}}^t\, (\,\,\,)\,d\tau $$ (eq. 12), case 2b traces the integral only for a fraction of the diffusion time, ds2/κs. For all cases, we compare our result with a formulation (Amiri & Vafai 1994; Minkowycz et al. 1999) in which the heat transfer between fluid and solid is described by the instantaneous temperature difference Tf–Ts and an interstitial heat transfer coefficient. Their formulation can be rewritten as   \begin{equation}{q_{{s_{{\rm{inst}}}}}}\! \left( {\vec{x},t} \right) = \frac{{\bar{\lambda }}}{{{\delta _M}}}\ \left( {{T_f}\left( {\vec{x},t} \right) - {T_s}\left( {\vec{x},t} \right)} \right),\end{equation} (27) where $$\bar{\lambda }$$ is the mean thermal conductivity of the solid and the fluid phases, in general an arithmetic ($$\bar{\lambda } = \ \varphi {\lambda _f} + ( {1 - \varphi } ){\lambda _s}$$) or geometric ($$\bar{\lambda } = \lambda _f^\varphi \ \cdot \lambda _s^{( {1 - \varphi } )}$$) average may be used, and δM is the microscopic thickness of the thermal boundary layer between the fluid and the solid phases. As indicated in eq. (27), the interstitial heat transfer coefficient can be understood as a mean thermal conductivity divided by an effective thermal boundary layer thickness δM. In the formulations by Minkowycz et al. (1999) and Amiri & Vafai (1994), the interstitial heat transfer coefficient is formulated in terms of constant model-dependent non-dimensional numbers such as the Sparrow number, Reynolds number and Prandtl number (for definitions, see their papers). From our formulation, however, it is clear that the effective thermal boundary layer is time dependent. It is thus interesting to test, to which extent constant values of δM may approximate the thermal non-equilibrium behaviour in porous flow. We use the subscript ‘inst’ for the resulting heat flow, since this formula gives the instant response to a temperature difference between the two phases, but does not consider the thermal history of the interphase heat exchange. The equations were programmed in MATLAB (MATLAB, 2012), using second-order finite differences (FDs) for the numerical equivalent of eqs (24) and (25), with upwind for the advection term in eq. (25), various integration rules for eq. (26) and a maximum Fourier coefficient nmax = 500 in eq. (23). Integration of eq. (26) needs special attention because the integrand has a singularity at the upper boundary due to the effective $${\sqrt {t'} ^{ - 1}}$$ behaviour of q*. In principle, for a finite number nmax of Fourier coefficients q* remains finite and the trapezoidal rule may be used to evaluate eq. (26), though with very small time steps (cf.Section 3.1). Alternatively, for case 2, we integrate eq. (26) using the Simpson rule up to the time t' − dt', where dt' is the time step (for an odd number of integration intervals we apply the Trapez rule for the very first interval, because the Simpson rule requires an even number of intervals). Rather than solving the remaining improper integral $$\mathop \smallint \nolimits_{t{'} - dt{'}}^{t{'}} A( {\tau {'}} )q_s^*( {t{'} - \tau {'}} )d\tau {'}$$ numerically we assume $$q_s^*( \tau )\ = \ a/\sqrt \tau $$ and utilize the analytic solution (using the mean value theorem) $$\mathop \smallint \nolimits_0^{d{t^{'}}} A( \tau )q_s^*( \tau )d\tau \ = \ 2A( {{\tau _{{\rm{mean}}}}} )q_s^*( {dt{'}} )dt{'}$$, that is, the argument A has to be taken at an intermediate time 0 ≤ τ'mean ≤ dt' and q*(dt') is taken at the first time increment dt'. In our calculations, we use τ'mean = 0.3dt' but other choices gave similar results. The numerical values of the non-dimensionalized variables are given in Table 2. Table 2. Numerical values of the non-dimensional variables. Symbol  Definition  Value  as  Geometry factor for the solid phase  2  $$d_s^{'}$$  Characteristic length of the solid phase  1  φ  Porosity  0.1–0.9  h'  Width of the model (depending on diffusion time)  0.125–10  $$T_{f0}^{'}$$  Initial temperature of the fluid phase  1  $$T_{s0}^{'}$$  Initial temperature of the solid phase  0  $${\rm{\Delta }}T_0^{'}$$  Initial Temperature difference between the fluid and the solid phases  1  δM'  Effective thermal boundary layer thickness  0.05–3  $$t_{{\rm{max}}}^{'}$$  Maximum time  0.5–3  dt'  Time step  0.1–0.00001  $$v_f^{'}$$  Velocity of the fluid (i.e. Peclet number)  0–300  dx'  Grid size of the (macroscopic) model  0.01–0.02  Symbol  Definition  Value  as  Geometry factor for the solid phase  2  $$d_s^{'}$$  Characteristic length of the solid phase  1  φ  Porosity  0.1–0.9  h'  Width of the model (depending on diffusion time)  0.125–10  $$T_{f0}^{'}$$  Initial temperature of the fluid phase  1  $$T_{s0}^{'}$$  Initial temperature of the solid phase  0  $${\rm{\Delta }}T_0^{'}$$  Initial Temperature difference between the fluid and the solid phases  1  δM'  Effective thermal boundary layer thickness  0.05–3  $$t_{{\rm{max}}}^{'}$$  Maximum time  0.5–3  dt'  Time step  0.1–0.00001  $$v_f^{'}$$  Velocity of the fluid (i.e. Peclet number)  0–300  dx'  Grid size of the (macroscopic) model  0.01–0.02  View Large 3.1 Case 1: solid and fluid phases in rest The first test case is an initial value problem with $$T_f^{'}( {t{'} = \ 0} ) = \ 1$$ and $$T_s^{'}( {t{'} = \ 0} ) = \ 0$$. Both the fluid and the solid are assumed at rest. On the macroscopic scale, no temperature variations are assumed, thus heat eqs (24) and (25) degenerate to   \begin{equation}\left( {1 - \varphi } \right)\ \frac{{\partial {T_s}{'}}}{{\partial t{'}}} = \ + {Q_b}^\prime \end{equation} (28) and   \begin{equation}\varphi \ \frac{{\partial {T_f}{'}}}{{\partial t{'}}} = \ - {Q_b}^\prime .\end{equation} (29) Before solving these equations with our method, we derive the analytical solution by switching to the microscopic scale. On the microscopic scale, the fluid is assumed to occupy infinite planar fluid layers, thus we have the rectangular temperature distribution discussed previously (see eq. 15) and Fig. 5). On the microscopic scale, the temperature in the fluid and the solid phases, respectively, can directly be evaluated by solving the non-dimensional form of eq. (15). To obtain the macroscopic fluid and solid temperatures, the obtained temperatures are averaged by spatial integration of eq. (15) over the fluid and solid volumes, respectively. The resulting temperatures represent the spatially constant fluid and solid temperatures on the macroscopic scale and are functions of time and of porosity (which is in our case the volumetric fraction of the fluid). The result is depicted in Fig. 7 and will be called in the following ‘the analytical solution’. Figure 7. View largeDownload slide Non-dimensional temperature versus non-dimensional time for the fluid phase (red) and the solid phase (blue) as a function of time for different porosities. The porosities φ valid for each curve are depicted at the end of the curves, outside of the box. The initial temperature disturbance is a rectangular half-wave. Figure 7. View largeDownload slide Non-dimensional temperature versus non-dimensional time for the fluid phase (red) and the solid phase (blue) as a function of time for different porosities. The porosities φ valid for each curve are depicted at the end of the curves, outside of the box. The initial temperature disturbance is a rectangular half-wave. For thin fluid layers (φ ∼ ≤ 0.3), the fluid is rapidly cooled, when coming into contact with the solid material, which in turn heats up until thermal equilibrium is reached at times t' ≤ 0.3. On the contrary, for wide fluid layers (φ ∼ ≥ 0.7), the temperature decrease in the fluid is much slower and equilibrium is reached after a much longer time. To test our approach, we calculate the IPHE term with the convolution integral (eq. 26), and solve the heat eqs (28) and (29) for porosities φ equal to 0.1, 0.5 and 0.9 and compare the results to the analytical solution and the approximate solutions based on eq. (27). The results are shown in Fig. 8. For each porosity case, we used three different non-dimensional time steps dt' for the integration of the convolution to show the convergence of our approach. The inset in Fig. 8(c) shows the convergence of the errors of the late equilibrated stages in a log–log diagram. For φ = 0.1 (Fig. 8a), the dykes are very thin and cool rapidly. On the microscopic scale, this behaviour is associated with a very steep temperature gradient in the boundary layer between the fluid and the solid phases. On the macroscopic scale, rapid cooling demands a sufficiently high temporal resolution, here small time steps of about 10−5 are needed, the error order is about 1/2 (Fig. 8a and inset in c). The approximate solution with instantaneous heat exchange (eq. 27) using δM = 0.35 underestimates the thermal equilibration in the fluid and the solid phases at early times (t' < 0.04), but gives good results for later times (Fig. 8d). For φ = 0.5 (Fig. 8b and inset in c), a time step of dt' = 10−4 gives a very good fit to the analytical solution with an error order also of about 1/2; the approximate solution still underestimates the thermal equilibration in the very early cooling stage, but is also very close to the analytical solution. Figure 8. View largeDownload slide Non-dimensional temperature for the solid phase (blue lines) and the fluid phase (red lines) versus non-dimensional time: (a)–(c) our approach: based on an initial rectangular half-wave temperature perturbation and (d)–(f) approximate solution: based on eq. (27), for porosities of (a) and (d) 0.1, (b) and (e) 0.5 and (c) and (f) 0.9 in comparison to the analytical solution (green lines, see Fig. 7). For (a)–(c) three different non-dimensional time steps dt΄ are used for each model to show the convergence of the approach; for (d)–(f) four different δ'M values are used (and the smallest time steps dt΄ according to the results in (a)-(c). Note the difference in t΄ scale. Figure 8. View largeDownload slide Non-dimensional temperature for the solid phase (blue lines) and the fluid phase (red lines) versus non-dimensional time: (a)–(c) our approach: based on an initial rectangular half-wave temperature perturbation and (d)–(f) approximate solution: based on eq. (27), for porosities of (a) and (d) 0.1, (b) and (e) 0.5 and (c) and (f) 0.9 in comparison to the analytical solution (green lines, see Fig. 7). For (a)–(c) three different non-dimensional time steps dt΄ are used for each model to show the convergence of the approach; for (d)–(f) four different δ'M values are used (and the smallest time steps dt΄ according to the results in (a)-(c). Note the difference in t΄ scale. For φ = 0.9 (Fig. 8c), the wide hot dykes heat up the thin solid phase in between the dykes. For this case, the convergence is poor (inset in Fig. 8c). Comparing Fig. 8(c) to (a), equilibration appears to be slow. However, since $$t{'} = \ t \cdot \frac{{\tilde{\kappa }}}{{d_s^2}}$$, t'scales with the quadratic inverse of the characteristic dykes distance ds which is different for the various φ − cases. In fact, the equilibration timescales with the diffusion time associated with the thinner layer, which, in non-dimensional form, is equal 1 in the case of φ = 0.9, but which is equal 0.1 in case of φ = 0.1. Using a time step of dt' = 10−3, we are very close to the analytical solution. Coarser time steps lead to an underestimation of the heat flow and thus, also underestimate the temperature in the solid for early times and lead to an overshot of the temperature at later times. The approximate solution (eq. 27) (Amiri & Vafai 1994; Minkowycz et al. 1999) gives good results for late times, but does not catch the thermal equilibration at early times, neither with $$\delta _M^{'}$$ = 0.35, nor with $$\delta _M^{'}$$ = 3.0. Fiddling with this parameter will certainly lead to an improved fit, and may allow the use of this approximative formula in certain cases, but the need to tune this parameter beforehand remains. 3.2 Case 2a: solid phase in rest, fluid phase moving In the next case, we assume that the fluid phase is moving with a constant velocity $$v_f^{'}$$ in z-direction. Due to non-dimensionalization, the velocity can directly be written as the Peclet number (i.e. $$v_f^{'} = \ {\rm{Pe}}$$) which can be defined by   \begin{equation}{\rm{Pe\ }} = \ \delta v \cdot {d_s} \cdot {\kappa ^{ - 1}},\end{equation} (30) with δv = vf − vs, i.e. in our case δv = vf. Since the z-direction is the only space coordinate, we can treat the problem in 1-D. We study the problem over a length scale of $$z_{{\rm{max}}}^{'} = \ 10$$. The initial conditions are $$T_f^{'}( {t{'} = \ 0} ) = T_s^{'}( {t{'} = \ 0} ) = \ 0$$ and the boundary conditions are $$T_f^{'}( {z{'} = \ 0} ) = T_s^{'}( {z{'} = \ 0} ) = \ 1$$ and $$T_f^{'}( {z{'} = z_{{\rm{max}}}^{'}} ) = T_s^{'}( {z{'} = z_{{\rm{max}}}^{'}} ) = \ 0$$. This condition implies a sudden heating at the bottom (z' = 0). The microscopic equivalent of this setup resembles Fig. 2(a) with hot fluid entering from below. Thus, our model describes the heating half-space problem including heat conduction and advection by porous flow. First, we explore the time evolution. Fig. 9(a) shows the solid and fluid temperatures assuming Pe = 3 and φ = 0.5 at times 0.25, 0.5, 1, 2 and 3. At the early times (0.25, 0.5), the pairs of curves resemble the erf function-like conductive cooling half-space solution, while at later times, advection takes over and the solutions develop into reverse s-shaped curves. This transitional behaviour stems from the fact that a conductively propagating temperature front moves with the velocity   \begin{equation}{v_{{\rm{diff}}}} = \ c\sqrt {\frac{\kappa }{t}} \end{equation} (31) with c of the order 1, while a convectively propagating front moves with the Darcy velocity $$\varphi v_f^{'}$$. Equating these velocities gives the overtake time, assuming c = 1)   \begin{equation}{t_{{\rm{ov}}}} = \frac{\kappa }{{v_f^2{\varphi ^2}}}\end{equation} (32) after which advection will take over and dominate the heat transport. This time can be written in non-dimensional form using the Peclet number to   \begin{equation}{t_{{\rm{ov}}}}^\prime = \frac{1}{{{\rm{P}}{{\rm{e}}^2}{\varphi ^2}}}\end{equation} (33) Figure 9. View largeDownload slide Temperature in the solid phase (blue lines) and the fluid phase (red lines) for successive times t'indicated in the figures, and for a setting with a porosity of 0.5. The fluid and the solid are suddenly heated up to T΄ = 1 at z΄ = 0 and the fluid moves upward in z’-direction with a velocity defined by the Peclet numbers (a) 3 and (b) 300. Figure 9. View largeDownload slide Temperature in the solid phase (blue lines) and the fluid phase (red lines) for successive times t'indicated in the figures, and for a setting with a porosity of 0.5. The fluid and the solid are suddenly heated up to T΄ = 1 at z΄ = 0 and the fluid moves upward in z’-direction with a velocity defined by the Peclet numbers (a) 3 and (b) 300. From this equation, it is clear that advection becomes important for the case shown in Fig. 9(a) for non-dimensional times larger than $${t_{{\rm{ov}}}}^{'} = \frac{4}{9}$$ (given Pe = 3 and φ = $$\frac{1}{2}$$), in agreement with the curve shapes changing from erf-type to S-type. Secondly, it is obvious, that from early times on the fluid is warmer than the solid. Thus, thermal non-equilibrium is clearly demonstrated. Hot fluid enters from below and is cooled by the solid, while the solid is heated by the advected fluid giving a temperature contrast along the curves of about 0.15 at early times decreasing slowly to 0.05 at t = 3. At later times, the temperature fronts of the fluid and solid move roughly with the same speed and the fluid temperature front is always about 0.25 diffusion lengths ahead. Thus, even on the macroscopic time and length scale such non-equilibrium temperatures are to be expected even for a Peclet number as low as 3. It is interesting to examine the evolution also in the limit of high Pe- numbers which is the relevant case for melt migrating through dykes. We carried out an experiment with Pe = 300 for φ = 0.5) and show the temperatures at the times 0.01, 0.02, 0.04, 0.08, 0.12, 0.16 and 0.2 in Fig. 9(b). One can clearly see the high-temperature front of the fluid entering and passing through the medium. Behind the front, the fluid has cooled only slightly to values above 0.9. The hot fluid acts as a continuous heating for the cold matrix within a region equivalent to the actual penetration depth of the fluid. Within this region, the matrix temperature continuously increases, while within the lowermost part (z' = 0–0.1) still a steep thermal boundary layer is present. Thus, a pronounced kink characterizes the resulting solid temperature profiles, separating the thin conductive thermal boundary layer and the deep region where penetrating fluid heats up the solid. Thus, for high Pe- numbers the maximum temperature difference between solid and fluid approaches 1. The lateral offset of the fluid and solid temperature curves increases indefinitely. In Fig. 10, we explore the effect of different Pe- numbers (fluid velocities) and porosities. The temperature in the fluid and the solid are shown at time t' = 1 (i.e. one diffusion time, see eq. 21) for an upward moving fluid for a porosity of φ = 0.1 (Fig. 10a) and φ = 0.5 (Fig. 10b), respectively. The results are shown for Pe- numbers equal to 0, 1, 3, 10 and 20 for φ = 0.5 (Fig. 10b) and for Pe- numbers equal 3, 10, 30 and 60 for φ = 0.1 (Fig. 10a). For the highest fluid velocity in each case (Pe = 60 for φ = 0.1 and Pe = 20 for φ = 0.5), the upper boundary condition at $$z_{{\rm{max}}}^{'}$$ strongly affects the results. We included the case for Pe = 0 (i.e. $$v_f^{'}$$ = 0), when the liquid is not moving, to verify that the solutions in the fluid and solid are identical and equal to the cooling (or heating) half-space solution (Turcotte & Schubert 2002) (as long as the heating front does not reach the right-hand side of the model). In fact, the analytical solution for the cooling half-space (which corresponds to Pe = 0) was also drawn in Fig. 10(b), and is perfectly matched by our solutions and thus not distinguishable. The time increments for the integration in eq. (26) depends on the porosity and were chosen according to the best solutions obtained for case 1 (see Fig. 8) to dt' = 10−5 for φ = 0.1 and to dt' = 10−4 for φ = 0.5. Figure 10. View largeDownload slide (a) and (b) Temperature in the fluid phase (red lines) and in the solid phase (blue lines) at the time t' = 1 (i.e. one diffusion time) for a setting where the fluid and the solid are suddenly heated up to T' = 1 at z' = 0 and the fluid is moving upward with different velocities (defined by the Peclet number Pe) for a porosity (a)φ = 0.1, (b)φ = 0.5, (c) non-dimensionalized temperature difference between fluid and solid phases and (d) spatial (i.e. vertical) offset of the fluid and solid temperature curves as a function of Peclet number times porosity (i.e. a measure for the volumetric flow rate). Note that the dimensionalized z is different in cases (a) and (b), because the scaling length ds is a function of φ. Figure 10. View largeDownload slide (a) and (b) Temperature in the fluid phase (red lines) and in the solid phase (blue lines) at the time t' = 1 (i.e. one diffusion time) for a setting where the fluid and the solid are suddenly heated up to T' = 1 at z' = 0 and the fluid is moving upward with different velocities (defined by the Peclet number Pe) for a porosity (a)φ = 0.1, (b)φ = 0.5, (c) non-dimensionalized temperature difference between fluid and solid phases and (d) spatial (i.e. vertical) offset of the fluid and solid temperature curves as a function of Peclet number times porosity (i.e. a measure for the volumetric flow rate). Note that the dimensionalized z is different in cases (a) and (b), because the scaling length ds is a function of φ. As can be expected, the temperature difference between the fluid and the solid phases increases with higher velocities (larger Peclet number) of the fluid. In case of a high volume fraction of the fluid (φ = 0.5, Fig. 10b), the heat content of the fluid is large and it heats the solid continuously up while moving upward. In case of a low volume fraction (φ = 0.1, Fig. 10a) its heat content is small and it is losing the heat quickly to the solid, rapidly cooling down, and the temperature in the large volume of the solid is only slightly increased. Thus higher Peclet numbers are needed to obtain a similar heating effect as for φ = 0.5. Interestingly, for similar penetration depths of the fluid and solid temperature fronts, their non-equilibrium temperature differences are similar and to first order independent of φ. The non-equilibrium temperature differences can be quantified by taking the temperature differences of the fluid and solid curves in Figs 10(a) and (b) at a given solid temperature (here 0.5) and plot the temperature difference (Fig. 10c) and the spatial distance between solid and fluid temperature values of 0.5 (Fig. 10d) as a function of the Darcy flow related Peclet number PeD  \begin{equation}{\rm{P}}{{\rm{e}}_D} = {\rm{Pe}} \cdot \varphi .\end{equation} (34) First, we observe that with increasing PeD the non-equilibrium temperature difference increases, but for the highest PeD the curves begin to flatten. It should be remembered that per definition the (non-dimensional) temperature difference is bounded by 1 and thus an indefinite increase is impossible (cf. Fig. 9b, which shows that for high Peclet number the temperature difference rapidly reaches values close to 1). Secondly, inspection of the time dependence of the temperature difference (e.g. shown in Fig. 9) indicates that the temperature difference decreases with time due to a decreasing temperature gradient with time. Thus, the curves shown in Fig. 10(c) show a snapshot, and for later times they will flatten. Thirdly, the spatial offset of the temperature curves (Fig. 10d) shows a clear linear relationship, depicted by a slope equal 1 in the log–log plot. Thus, for large PeD the spatial offset increases unboundedly. This is particularly important, as from inspection of the time behaviour (Fig. 9a), the offset of the curves does not show a detectable time-dependences, indicating that (i) even for times large compared to one diffusion time the spatial offset is still present, and (ii) the offset is approximately independent of the overall temperature gradient. 3.3 Comparison to an FE solution We also tested our formulation against a solution obtained from an explicit model of a propagating hot fluid through a set of parallel sheets in a colder surrounding. For this purpose, we created a finite-element model (FE, using COMSOL) representing a vertical section of non-dimensional length 10 of a fluid channel embedded in a resting solid. This setup represents case 2a on the local scale. Only half of one fluid layer in contact with half of a surrounding solid layer are modeled (with a resolution of 50 × 500), assuming symmetric thermal boundary conditions at both sides. The initial temperature was a set to zero everywhere. At the initial time, a temperature jump of amplitude 1 was prescribed at the lower boundary to both solid and fluid boundary sections. The fluid was assumed to move upwards with a constant non-dimensional velocity of Pe = 20 within the fluid layer. At the upper boundary, we used a zero-temperature gradient boundary condition both for the solid and fluid boundary sections. The thermal conductivities of fluid and solid are assumed to be equal. For this system, we solved the energy conservation equation, $$\frac{{\partial T{'}}}{{\partial t{'}}} + {\vec{v}_f}{'} \cdot \vec{\nabla }T{'}\ = \ {\nabla ^2}T{'}$$, with $${\vec{v}_f}{'}$$ abruptly changing to zero at the fluid–solid interface. The fluid fraction (i.e. porosity) was 50 per cent. The resulting temperature fields within the fluid and the solid phases are shown in Fig. 11(a) for different times. The hot fluid penetrates into the cold solid and generates a wavy- or finger-like hot temperature front which is diffusively broadening vertically with time, while the amplitude of the wave remains roughly constant. This is equivalent to the approximate time independence of the spatial distances of the fluid–solid temperature curves, for example, in Fig. 9(a). Figs 11(b)–(d) show the horizontally averaged temperatures within the solid and fluid phases in comparison with the temperatures determined by the macroscopic approach (eqs 24–26) using full integration. The agreement is very good. At the temperature front the numerical solutions seem to exhibit little numerical diffusion (see also Fig. 11a). At all times, the temperature differences between the numerical and macroscopic solutions are small, always less than 2 per cent. These differences result from both the resolution errors of the numerical FE model and the FD approximation of the heat eqs (23) and (25), as well as the numerical integration in eq. (26). Figure 11. View largeDownload slide (a) Three snapshots of the non-dimensional temperature field of an FE solution on the local scale. The computation domains have been mirrored four times for better visualization. (b)–(d) Comparison between the horizontally averaged temperatures (fluid and solid) obtained by the macroscopic approach as described in Section3.2 and by the FE model for a fluid sheet in a solid surrounding, resolving the fluid and the solid explicitly. Figure 11. View largeDownload slide (a) Three snapshots of the non-dimensional temperature field of an FE solution on the local scale. The computation domains have been mirrored four times for better visualization. (b)–(d) Comparison between the horizontally averaged temperatures (fluid and solid) obtained by the macroscopic approach as described in Section3.2 and by the FE model for a fluid sheet in a solid surrounding, resolving the fluid and the solid explicitly. 3.4 Comparison to approximate formula We also calculated the temperatures in both phases for a range of Peclet numbers (between 3 and 100) and fluid fractions φ (0.1 and 0.5) based on the approximate formula (eq. 27) and for different values of $$\delta _M^{'}$$ (i.e. the thickness of the microscopic thermal boundary layer) and compared the results to those obtained by the full backward in time integration. Fig. 12 shows the cases for (Pe, φ) equal to (30, 0.1) and (10, 0.5) and $$\delta _M^{'}$$ values of 0.05, 0.2, 0.3 and 0.5. Inspecting Fig. 12 suggests for the first case a $$\delta _M^{'}$$ value of about 0.2 and in the second case a value around 0.3. We calculated the rms deviations between the full integral solution and the approximate solutions for all cases and determined the minimum for the rms curves, which was $$\delta _M^{'}$$ (Pe = 30, φ = 0.1) = 0.2102 for the first case and $$\delta _M^{'}$$ (Pe = 10, φ = 0.5) = 0.3352 for the second case. As can be expected, in general we found a decrease in $$\delta _M^{'}$$ with higher Peclet numbers (in agreement with Amiri & Vafai 1994). However, we also found cases for small Peclet numbers where $$\delta _M^{'}$$ did not depend on Pe. Figure 12. View largeDownload slide Temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 based on the approximate eq. (27) (thick lines) and different values of $$\delta _M^{'}$$ in comparison with the results of the full backward integration in time (thin lines; same results as in Figs 10a and b) for fluid fraction φ = 0.1 and Peclet number Pe = 30 (upper row) and φ = 0.5 and Pe = 10 (lower row). Figure 12. View largeDownload slide Temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 based on the approximate eq. (27) (thick lines) and different values of $$\delta _M^{'}$$ in comparison with the results of the full backward integration in time (thin lines; same results as in Figs 10a and b) for fluid fraction φ = 0.1 and Peclet number Pe = 30 (upper row) and φ = 0.5 and Pe = 10 (lower row). To test whether the optimum $$\delta _M^{'}$$ changes with time (cf. Fig. 9 showing the time evolution of the full solution), we determined the optimum $$\delta _M^{'}$$ (Pe = 10, φ = 0.5) as 0.3352 at time t' = 1 and ranan approximate case with this $$\delta _M^{'}$$ for the time 0–3. The resulting evolution agrees extremely well (within line width) with the full solution (Fig. 9) for all times examined. This might be surprising because in case 1 we observed a significant disagreement between the approximate and the exact solution, however, this disagreement is most pronounced only for early times (t' < 0.25), while for later times as examined in Fig. 9, the instantaneous non-equilibrium approach seems to be sufficient if an appropriate $$\delta _M^{'}$$ is chosen. To summarize, the approximate formula (eq. 27) is a good approximation (and much simpler to implement in a numerical code for melt segregation) with a $$\delta _M^{'}$$value between 0.2 and 0.3 (depending on the average fluid (melt) velocity). 3.5 Case 2b: testing the necessary length of the convolution integral So far, we have carried out the integration in eq. (26) from the initial time $$t_0^{'}$$ = 0 to t' (i.e. along the entire fluid path). Inspecting the shape of the heat flow response function q* (Fig. 6) which rapidly decays for times of one diffusion time or less it might be expected that it is not necessary to carry out the convolution integral for the full time of evolution in order to save computation time. To check this, we defined a certain truncation time $$t_{{\rm{trunc}}}^{'} = t{'}\ - \Delta t_{{\rm{trunc}}}^{'}$$ and tested the quality of the solutions for the case φ = 0.1, a Peclet number of 3 and 10, and $$\Delta t_{{\rm{trunc}}}^{'}$$ = 0., 10−5, 10−4, 10−3, 10−2, 10−1 and 1. Results are plotted in Fig. 13 and show that the temperature in the fluid is strongly overestimated if the backward-in-time integration in eq. (26) is too short, while the temperature in the solid is much less affected. We conclude that always the full backward in time integrations is needed to avoid errors in solid–fluid heat transfer. The application to generalized two-phase flow requires a Lagrangian formulation to evaluate the integral and is thus computationally expensive. However, to avoid repeated calculations, it is advisable to store the result of the parentheses in the integrals in eqs (11) or (12) during the forward stepping. Figure 13. View largeDownload slide Dependency of the accuracy of the result on the length of the backward-in-time integration in eq. (26) for (a) slow fluid velocity (Pe = 3) and (b) faster fluid velocity (Pe = 10). Shown is the temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 and for a porosity of φ = 0.1. Figure 13. View largeDownload slide Dependency of the accuracy of the result on the length of the backward-in-time integration in eq. (26) for (a) slow fluid velocity (Pe = 3) and (b) faster fluid velocity (Pe = 10). Shown is the temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 and for a porosity of φ = 0.1. 4 DISCUSSION AND CONCLUSION 4.1 General remarks In this paper, we developed a theory to determine the heat exchange between a fluid and a solid phase on a macroscopic scale where the fluid is intruding or segregating into the solid and fluid and solid have initially different temperatures. The theory allows us to test the conditions under which different temperatures may prevail in the fluid and the solid phases. Such a formulation is a pre-requisite to model melt movement within solid rocks in regions where the temperature of the solid rock is below the solidus temperature. It may also be important when describing hydrothermal groundwater flow and heat transport through large fractures, which are not known well enough to be resolved as single structures. Heat exchange is based on a formula which implicitly assumes fluid-filled layers of equal spacing where the width of the layers depends on the volumetric fluid fraction. For a moving solid phase with a differentially moving fluid phase we derived an expression to calculate the heat exchange rate across the fluid–solid interfaces using this formula and integrating the heat exchange history along solid flow path. Solving this integral gives an interphase heat exchange rate of dimension Jm−3 s−1, which couples the two heat transport equations for the solid and the fluid phases. Solving these heat equations allows determining the spatially and temporarily varying temperatures in the fluid and the solid phases separately. 4.2 Testing thermal non-equilibrium by using Peclet numbers Before discussing the possible applications of our formulation, it is worthwhile to introduce alternative Peclet numbers which may be more appropriate for various scenarios. While our Peclet number (eq. 30) is based on the fluid velocity and a length scale describing the distance ds of fluid inclusions (dykes and pores), a macroscopic Peclet number based on the length scale of the system L may be defined by   \begin{equation}{\rm{P}}{{\rm{e}}_L} = \frac{{\delta vL}}{\kappa }\quad {\rm{with}}\quad {\rm P}{{\rm e}_L} = \ {\rm{Pe}}\frac{L}{{{d_s}}}\end{equation} (35) Often the volumetric flow rate (Darcy velocity vD) is used instead of the fluid velocity leading to the Darcy related Peclet number   \begin{equation}{\rm{P}}{{\rm{e}}_D} = \frac{{{v_D}{d_s}}}{\kappa }\quad {\rm{with}}\quad {\rm{P}}{{\rm{e}}_D} = \ {\rm{Pe}}\ \varphi \quad {\rm{from}}\,{\rm{eq}}.\,\,(34).\end{equation} (36) It is noteworthy that the overtaking time tov' at which Darcy related advection takes over with respect to heat conduction (cf. eq. 33) is directly given by PeD− 2. If the Peclet number is based on the Darcy velocity and a macroscopic length scale of the system L, we obtain   \begin{equation}{\rm{P}}{{\rm{e}}_{LD}} = \frac{{{v_D}L}}{\kappa }\quad {\rm{with}}\quad {\rm{P}}{{\rm{e}}_{LD}} = \ {\rm{Pe}}\ \frac{{\varphi L}}{{{d_s}}}\end{equation} (37) In many applications of two-phase flow, neither the microscopic, nor the macroscopic length scales are used, but the problems are scaled by the compaction length (McKenzie 1984). Using McKenzie's definition, the compaction length δc can be written in terms of our characteristic length scale ds by   \begin{equation} {\delta _c} = {\left( {\frac{{\left( {{\eta _b} + \frac{4}{3}{\eta _s}} \right){\varphi ^n}}}{{{\eta _f}b}}} \right)^{{\raise0.7ex\hbox{1} / \lower0.7ex \hbox{2}}}}\ {d_s}. \end{equation} (38) Here ηb, ηs and ηf are the bulk, shear and fluid viscosities, respectively, b is a geometric constant of the order 100 and n is the power of the permeability–porosity relation. For realistic partially molten systems, the compaction length is thus 2–7 orders of magnitude larger than ds. The main results of this paper may be used to check for a given Darcy flow scenario whether the temperature field is dominated by Darcy advection and whether fluid and solid temperatures are different due to thermal non-equilibrium. A quick check starts from determining the characteristic Peclet number and porosity of the system to obtain PeD. If the characteristic non-dimensional timescale of the flow problem is larger than PeD− 2, then the temperature field will be dominated by Darcy flow rather than conduction alone. Furthermore, Figs 10(c) and (d) imply that in the presence of an overall temperature gradient in fluid flow direction fluid and solid are not in thermal equilibrium if PeD > O(1) to O(10). Under such conditions, the temperatures may be different by 10’s per cent of the typical temperature variations of the problem. This is equivalent to the fact that the same temperatures in the solid and fluid phases may be spatially separated by a distance which is proportional to the characteristic length of the solid phase, ds, times PeD. 4.3 Combining the thermal non-equilibrium formulation with general two-phase flow In the test cases as presented in this paper, the fluid velocity is prescribed. However, in a hydrothermal flow problem with a given rock permeability distribution, the fluid flow can be calculated by solving additionally the Darcy and mass conservation equations. The convolution integral over time (eq. 26) still needs to be evaluated for every point in space. In the full two-phase flow case, the matrix may be deformable and may flow according to solutions of the Stokes equations. To apply our thermal non-equilibrium formulation first, the flow paths of the matrix (not of the fluid) need to be determined and stored. At a given time step, the convolution integral must be carried out backward in time along the matrix path, which implies that the temperature changes need to be known everywhere for all times before the given time step. While rigorous thermal non-equilibrium modeling of porous flow requires the numerically expensive full backwards integrations of the convolution integral (which means that the temperature fields for solid and fluid phases have to be stored for each elapsed time step and considered in the convolution integral for each new time step), we showed that for times larger than the characteristic diffusion time the approximate formula (eq. 27) (cf. approximations by Amiri & Vafai 1994 or Minkowycz et al. 1999) does a good job and significantly reduces computation times. This is particularly important for full two-phase flow formulations in which our theory would imply that also the evolving flow paths of the matrix need to be tracked and integrated backward in time. Thus, if this integration could be replaced by an approximate instantaneous heat exchange term, two-phase flow formulations such as shown by, for example, Schmeling (2000) or Wallner and Schmeling (2016) could be easily extended by separating the heat equation into two equations and coupling them by the instantaneous heat exchange term. Our method may also be applied to fluid geometries others than regularly spaced channels or dykes. This requires alternative formulations of the heat flow response function to a temperature step function, $$q_s^*( t ),$$ that is, replacement of eq. (17). For more complicated geometries, such functions can be determined numerically by solving an initial value problem of heat conduction and storing $$q_s^*( t )$$ in a lookup table. We also want to emphasize that thermal non-equilibrium in hydrothermal porous flow problems is generally of transient nature, arising from a sudden change in the physical state: for example, a sudden change in the flow as when a dyke is injected. Thus, our formulation needs to be lined to a formalism which calculates of estimates the changing flow field. In case of steady state or slowly changing flow fields, sudden changes in the physical state may include porous-driven temperature fronts, in which thermal non-equilibrium is predicted to occur in the high-temperature gradient regions as shown in this paper. 4.4 Application to melt segregation at the transition between channeling and dykes One of the most obvious two-phase flow systems where our theory may be applied is the transitional regime between melt channeling within a melting source zone (Stevenson 1989; Aharonov et al. 1995) and magma ascent by dykes through the subsolidus mantle lithosphere or crust (e.g. Rivalta et al.2015). Keller et al. (2013) studied isothermal melt ascent by numerical models, focusing on rheology but neglected freezing, and identified regimes of viscous diapirism, viscoplastic decompaction channels and elastoplastic dyking. Leitch and Weinberg (2002) modeled mesoscale pervasive magma ascent through dykes but assumed thermal equilibrium between dykes and country rocks. Our results allow testing this kind of assumption. We will make a short, incomplete tour through these processes with the emphasis of extracting the parameters controlling the Darcy related Peclet number PeD, that is, we need estimates about the channel or dyke distance, ds and the volumetric flow rate vD = vf φ. See also Fig. 14 for illustration of this tour. Figure 14. View largeDownload slide Illustration of the discussion of melt ascent and associated Peclet numbers which result in predictions about thermal equilibrium or non-equilibrium during this ascent process. The horizontal solidus temperature indicates the situation, where melt temperatures will be slightly below solidus temperatures, thus the melt has the potential to penetrate into the subsolidus region. Figure 14. View largeDownload slide Illustration of the discussion of melt ascent and associated Peclet numbers which result in predictions about thermal equilibrium or non-equilibrium during this ascent process. The horizontal solidus temperature indicates the situation, where melt temperatures will be slightly below solidus temperatures, thus the melt has the potential to penetrate into the subsolidus region. From the analysis of Stevenson (1989), a partially molten rock undergoing deviatoric straining develops a channeling instability with a characteristic wavelength equal or smaller (down to several grain sizes) than the compaction length, which is defined by $${( {{k_\varphi }( {{\eta _b} + \frac{{4{\eta _s}}}{3}} )/{\eta _f}} )^{1/2}}$$, with kφ as initial permeability, ηb,  and  ηs as bulk and shear viscosity of the matrix and ηf as melt viscosity. Typical values for partially molten mantle give compaction lengths between 10s of metres and few kilometres. Models on channel formation by Richardson (1998) and Golabek et al. (2008) found channeling instabilities on larger scales, but these were probably limited to the finite grid resolution, inhibiting the development of smaller wavelengths. Well-resolved models of vein formation due to shearing and compaction by Rabinowicz & Toplis (2009) found the formation of regular veins with distances of 0.03 compaction lengths (several 10s of metres). Modeling of channeling instabilities due to reaction infiltration in compacting media by Spiegelman et al. (2001) found initial distances between veins of about 0.3 compaction lengths, but upon ascending they continuously merge into fewer channels with larger distances. Such coalescence of melt conduits has been proposed as the ‘fractal tree model’ by Hart (1993) based on geochemical arguments. Observation (e.g. Weinberg et al. 2015) shows that the distance between veins may be on the order of decimetres. In summary, we expect that the channel distance close to the onset of channeling within the melt source regions may be rather small, on the order of decimetres, but further up, still in the supersolidus partially molten regions, may increase up to the metre or kilometre scale. The segregation velocity of magma in highly permeable or fully melt-filled veins depends mainly on the melt viscosity. A velocity of around 2.5 m yr−1 (Petford 1995) has been estimated as an average value, however, this will strongly depend on the spacing and widths of the channels. We might therefore better try to estimate the volumetric flow rates (Darcy velocity). Let us take as a typical example of melt ascent a mid-ocean ridge, where a basaltic crust of thickness h spreads with a half-spreading rate of v0. If the melt is collected within an asthenospheric wedge of half-width L (which, as a consequence of the wedge shape, decreases towards the surface), on first order the Darcy velocity within the mantle wedge is given by v0h/L. For typical mid-ocean ridges, this results in a volumetric melt flow rate of 1e-10 to 1e-9 ms−1, the smaller value being more typical for greater depth (where L is wider). Such values have also been reported from the fully dynamic melting models by Keller et al. (2017). With a thermal diffusivity κ = 10−6 m2 s−1, a small (because deep) Darcy velocity and dense veins, the resulting Pe- numbers (PeD) are of the order of 10−5–10−3, that is, the flow is in thermal equilibrium. At shallower depth, and when melt conduits have merged, PeD may increase to order 1 or larger, that is, in the presence of a vertical temperature gradient thermal non-equilibrium effects will become moderately important (see Fig. 14). It should be noted that these rough estimates are conservative estimates, because they assume (1) steady-state ascent of melt and (2) passive upwelling. Regions with active mantle upwelling and time-dependent melt extraction will exhibit higher Peclet numbers implying a high potential for thermal non-equilibrium. When the above-mentioned melt channels reach the level of the solidus temperature, they must undergo a transition into dykes. Dykes are generally believed to be generated by magmatic overpressure and brittle failure of rocks and propagate through initially unfractured rock (Rivalta et al. 2015). We will now focus on this end-member of magma ascent and estimate appropriate Peclet numbers. The sizes of dykes in nature is very different, however, a typical dyke has a thickness between some 10s of centimetres to some 10s of metres with a length of a few 100 m. In the vertical plane, dykes (or in the horizontal plane for sills) can be some tens of kilometres (Wright et al. 2006). There is also a large amount of literature reporting the occurrence of tens to thousands of subparallel dykes and/or melt veins. Such dykes can be found, for example, in the vicinity of active volcanoes (Gudmundsson 1983; 2002), in ophiolites formed near mid-oceanic ridges (e.g. Nicolas et al. 1994) or in continental volcanic provinces (e.g. Ray et al. 2007) and the magmatic sources are proposed in the upper or lower crust (e.g. Guernina & Sawyer 2003) or even in the upper mantle (Ernst et al. 2001). Whether multiple dykes in a swarm are propagating at the same time is still a matter of investigation (e.g. Jin & Johnson 2008). Havlin et al. (2013) modeled dyke initiation from a melt accumulation layer at the top of the asthenosphere and obtained dyke recurrence times of the order of days and spacings of 100 m or larger. Field observations show at least that magma propagation in adjacent dykes during active rifting events occur within short time intervals of a few months (e.g. Einarsson & Brandsdottir 1980; Grandin et al. 2010). Dahm (2000) estimated for an average dyke with a width of a few decimetre and a length of a few kilometres a propagation velocity (which is equivalent with the propagation of the hot magmatic fluid) on the order of 0.1–1 ms−1: recalling the definition of the Peclet number (eq. 30) and assuming a dyke spacing of 100 m leads to a Peclet number on the order of 105, thus by far larger than studied in this papers, indicating a strong non-equilibrium temperature distribution between fluid and solid. From this exercise, we conclude (Fig. 14) that within the source regions of mantle melting steady-state porous flow through a channelized system is in thermal equilibrium with the matrix. Once this flow becomes episodic, with volumetric flow rates exceeding the steady-state rate be one or two orders of magnitude, the Peclet number becomes so large that thermal non-equilibrium effects may become important. Furthermore, within the transition zone between melt rising through highly permeable channels and melt rising through propagating dykes the Peclet number is expected to strongly exceed the value one. Within this transition zone, our results indicate that melt and matrix temperature may be different by an order of 10 per cent of the characteristic temperature variations of the tectonic setup, implying that melt is still above the solidus, while the solid is below. Entering the dyke regime, extremely large Peclet numbers are predicted implying large temperature differences between rock and fluid. 4.5 Summary We conclude with summarizing the main key points of this study: A formulation of calculating thermally non-equilibrium two-phase flow has been developed based on the time-dependent heat exchange between fluid and solid which integrates the thermal history along flow paths by a convolution integral. For Darcy velocity-based Peclet numbers larger 1 thermal non-equilibrium between fluid and matrix becomes important with temperature differences of 5 per cent or more relative to characteristic temperature differences of the flow problem. The numerically expensive calculation of the history-dependent heat exchange can be replaced in the long-term limit by using an approximate instantaneous exchange term with adequate parameters. Our results imply that thermal non-equilibrium can play an important role for melt migration through partially molten systems where melt focuses into melt channels and near the transition to melt ascent by dykes. Acknowledgements This paper was accomplished during a sabbatical stay of the authors at the Monash University, Melbourne. We are grateful to all our colleagues from the School of Earth, Atmosphere and Environment for their hospitality and support. We thank two anonymous reviewers for their patience to go through all equations and to contribute several constructive comments, which helped us to improve the paper considerably. We acknowledge funding support for this sabbatical stay by the Deutsche Forschungsgemeinschaft (DFG) with the grant no. SCHM 872/25-1. REFERENCES Aharonov E., Whitehead J.A., Kelemen P.B., Spiegelman M., 1995. Channeling instability of upwelling melt in the mantle. J. geophys. Res. , 100( 20), 433−450. Amiri A., Vafai K., 1994. Analysis of dispersion effects and non-thermal equilibrium, non-Darcian, variable porosity incompressible flow through porous media, Int. J. Heat Mass Transf. , 37 939– 954. Google Scholar CrossRef Search ADS   Bercovici D., Ricard Y., Schubert G., 2001. A two-phase model for compaction and damage, Part 1: general theory, J. geophys. Res. , 106( B5), 8887– 8906. Google Scholar CrossRef Search ADS   Dahm T., 2000. On the shape and velocity of fluid-filled fractures in the Earth, Geophys. J. Int . 142 181– 192. Google Scholar CrossRef Search ADS   de Lemos M.J.S., 2016. Thermal Non-equilibrium in Heterogeneous Media , Springer Science+Business Media, Inc. Google Scholar CrossRef Search ADS   Golabek G.J., Schmeling H., Tackley P.J., 2008. Earth's core formation aided by flow channeling instabilities induced by iron diapirs, Earth planet. Sci. Lett.   271 24– 33. Google Scholar CrossRef Search ADS   Grandin R., Socquet A., Jacques E., Mazzoni N., de Chebalier J.B., King G., 2010. Sequence of rifting in Afar, Manda-Hararo rift, Ethiopia, 2005–2009: time–space evolution and interactions between dikes from interferometric synthetic aperture radar and static stress change modeling, J. geophys. Res. , 115 B10413, doi:10.1029/2009JB000815. Google Scholar CrossRef Search ADS   Gudmundsson A., 1983. Form and dimensions of dykes in eastern Iceland, Tectonophysics , 95 295– 307. Google Scholar CrossRef Search ADS   Gudmundsson A., 2002. Emplacement and arrest of sheets and dykes in central volcanoes, J. Volc. Geotherm. Res. , 116 279– 298. Google Scholar CrossRef Search ADS   Guernina S., Sawyer E.W., 2003. Large-scale melt-depletion in granulite terranes: an example from the Archean Ashuanipi Subprovinces of Quebec, J. Metamorphic Geol. , 21 181– 201. Google Scholar CrossRef Search ADS   Einarsson P., Brandsdottir B. 1980. Seismological evidence for lateral magma intrusion during the July 1978 deflation of the Krafla volcano in NE Iceland, J. Geophys. , 47 160– 165. Ernst R.E., Grosfils E.B., Mège D., 2001. Giant dike swarms: Earth, Venus, and Mars. Annu. Rev. Earth planet. Sci. , 29 489– 534. Google Scholar CrossRef Search ADS   Hart S.R., 1993. Equilibrium during melting: a fractal tree model, Proc. Natl. Acad. Sci. USA , 90 11914– 11918. Google Scholar CrossRef Search ADS   Havlin C., Parmentier E.M., Hirth G., 2013. Dike propagation driven by melt accumulation at the lithosphere asthenosphere boundary, Earth planet. Sci. Lett. , 376 20– 28. Google Scholar CrossRef Search ADS   Jin Z.-H., Johnson S.E.. 2008. Magma-driven multiple dike propagation and fracture toughness of crustal rocks, J. geophys. Res. , 113 B03206, doi:10.1029/2006JB004761. Google Scholar CrossRef Search ADS   Katz R.F., 2008. Magma dynamics with the enthalpy method: benchmark solutions and magmatic focusing at mid‐ocean ridges, J. Petrol. , 49( 12), 2099– 2121. Google Scholar CrossRef Search ADS   Keller T., May D.A., Kaus B.J.P., 2013. Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust, Geophys. J. Int. , 195( 3), 1406– 1442. Google Scholar CrossRef Search ADS   Keller T., Katz R.F., Hirschmann M.M., 2017. Volatiles beneath mid-ocean ridges: deep melting, channelised transport, focusing, and metasomatism, Earth planet. Sci. Lett. , 464 55– 68. Google Scholar CrossRef Search ADS   Leitch A.M., Weinberg R.F., 2002. Modellinggranite migration by mesoscale pervasive flow, Earth planet. Sci. Lett. , 200 131– 146. Google Scholar CrossRef Search ADS   MATLAB and Statistics Toolbox Release 2012b, The MathWorks, Inc., Natick, MA. McKenzie D., 1984. The generation and compaction of partially molten rock, J. Petrol. , 25 713– 765. Google Scholar CrossRef Search ADS   Minkowycz W.J., Haji-Sheikh A., Vafai K., 1999. On departure from local thermal equilibrium in porous media due to a rapidly changing heat source: the Sparrow number, Int. J. Heat Mass Transf. , 42 ( 18), 3373– 3385. Google Scholar CrossRef Search ADS   Nield D.A., Bejan A., 2006. Convection in Porous Media , 3rd edn, Springer Science+Business Media, Inc. Nicolas A., Boudier F., Ildefonse B., 1994. Evidence from the Oman ophiolite for active mantle upwelling beneath a fast spreading ridge, Nature , 370 51– 53. Google Scholar CrossRef Search ADS   Petford N., 1995. Segregation of tonalitic-trondhjemitic melts in the continental crust: the mantle connection, J. geophys. Res. , 100 15735– 15743. Google Scholar CrossRef Search ADS   Rabinowicz M., Toplis M.J., 2009. Melt segregation in the lower part of the partially molten mantle zone beneath an Oceanic Spreading Centre: numerical modelling of the combined effects of shear segregation and compaction, J. Petrol. , 50 ( 6), 1071– 1106. Google Scholar CrossRef Search ADS   Ray R., Sheth H.C., Mallik J., 2007. Structure and emplacement of the Nandurbar–Dhule mafic dyke swarm, Deccan Traps, and the tectonomagmatic evolution of flood basalts, Bull. Volc ., 69 537– 551. Google Scholar CrossRef Search ADS   Richardson C.N., 1998. Melt flow in a variable viscosity matrix, Geophys. Res. Lett. , 25 1099– 1102. Google Scholar CrossRef Search ADS   Rivalta E., Taisne B., Buger A.P., Katz R.F. 2015. A review of mechanical models of dike propagation: schools of thought, results and future directions, Tectonophysics , 638 1– 42. Google Scholar CrossRef Search ADS   Rubin A.M., 1998. Dike ascent in partially molten rock, J. geophys. Res. , 103 ( B9), 20901– 20919. Google Scholar CrossRef Search ADS   Rudge J.F., Bercovici D., Spiegelman M., 2011. Disequilibrium melting of a two phase multicomponent mantle, Geophys. J. Int. , 184( 2), 699– 718. Google Scholar CrossRef Search ADS   Schmeling H., 2000. Partial melting and melt segregation in a convecting mantle, in Physics and Chemistry of Partially Molten Rocks , eds. Bagdassarov N., Laporte D., Thompson A.B., pp. 141– 178. Kluwer Academic Publ., Dordrecht. Google Scholar CrossRef Search ADS   Spiegelman M., 1993. Flow in deformable porous media. Part 1. Simple analysis, J. Fluid Mech. , 247 17– 38. Google Scholar CrossRef Search ADS   Spiegelmann M., Kelemen P.B., Aharonov, E., 2001. Causes and consequences of flow organization during melt transport: the reaction infiltration instability in compactible media, J. geophys. Res. , 106 2061– 2077. Google Scholar CrossRef Search ADS   Stevenson D.J., 1989. Spontaneous small-scale melt segregation in partial melts undergoing deformation, Geophys. Res. Lett. , 16( 9), 1067– 1070. Google Scholar CrossRef Search ADS   Turcotte D.L., Schubert G., 2002. Geodynamics . Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Wallner H., Schmeling H., 2016. Sensitivity analysis of rift induced delamination with application to Rwenzori Mountains, Geophys. J. Int. , 187 1135– 1145. Google Scholar CrossRef Search ADS   Weinberg R.F., Veveakis E., Regenauer-Lieb K., 2015. Compaction-driven melt segregation in migmatites, Geology , doi:10.1130/G36562.1. Wright T., Ebinger C., Biggs J., Ayele A., Yirgu G., 2006. Magma-maintained rift segmentation at continental rupture in the 2005 Afar dyking episode, Nature , 442 291– 294. Google Scholar CrossRef Search ADS PubMed  © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

A porous flow approach to model thermal non-equilibrium applicable to melt migration

Loading next page...
 
/lp/ou_press/a-porous-flow-approach-to-model-thermal-non-equilibrium-applicable-to-ADrRTpeL7Z
Publisher
Oxford University Press
Copyright
© The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx406
Publisher site
See Article on Publisher Site

Abstract

Abstract We develop an approach for heat exchange between a fluid and a solid phase of a porous medium where the temperatures of the fluid and matrix are not in thermal equilibrium. The formulation considers moving of the fluid within a resting or deforming porous matrix in an Eulerian coordinate system. The approach can be applied, for example, to partially molten systems or to brine transport in porous rocks. We start from an existing theory for heat exchange where the energy conservation equations for the fluid and the solid phases are separated and coupled by a heat exchange term. This term is extended to account for the full history of heat exchange. It depends on the microscopic geometry of the fluid phase. For the case of solid containing hot, fluid-filled channels, we derive an expression based on a time-dependent Fourier approach for periodic half-waves. On the macroscopic scale, the temporal evolution of the heat exchange leads to a convolution integral along the flow path of the solid, which simplifies considerably in case of a resting matrix. The evolution of the temperature in both phases with time is derived by inserting the heat exchange term into the energy equations. We explore the effects of thermal non-equilibrium between fluid and solid by considering simple cases with sudden temperature differences between fluid and solid as initial or boundary conditions, and by varying the fluid velocity with respect to the resting porous solid. Our results agree well with an analytical solution for non-moving fluid and solid. The temperature difference between solid and fluid depends on the Peclet number based on the Darcy velocity. For Peclet numbers larger than 1, the temperature difference after one diffusion time reaches 5 per cent of $$\tilde{T}$$ or more ($$\tilde{T}$$ is a scaling temperature, e.g. the initial temperature difference). Thus, our results imply that thermal non-equilibrium can play an important role for melt migration through partially molten systems where melt focuses into melt channels near the transition to melt ascent by dykes. Our method is based on solving the convolution integration for the heat exchange over the full flow history, which is numerically expensive. We tested to replace the heat exchange term by an instantaneous, approximate term. We found considerable errors on the short timescale, but a good agreement on the long timescale if appropriate parameters for the approximate terms are used. We derived these parameters which may be implemented in fully dynamical two-phase flow formulations of melt migration in the Earth. Hydrothermal systems, Mantle processes, Heat generation and transport, Mid-ocean ridge processes, Magma genesis and partial melting 1 INTRODUCTION Melt segregation from the deeper Earth to the upper crust or hydrothermal fluid flow in sedimentary basins are important examples of two-phase flows in geoscientific studies. Numerical modeling of such settings to quantify fluid flow and heat transport are commonly based on the Darcy equation describing fluid flow through a porous medium on a macroscopic scale (i.e. flow paths are not resolved). Porosity/permeability may be an intrinsic property of the rock or might be created by the formation of melt. Fluid flow is driven by thermal or compositional density differences between the fluid and the solid phases and a pressure gradient related to the state of stress of the fluid-filled rock. This typically applies to melt segregation within a partially molten rock volume, for example, beneath a mid-ocean ridge or within a plume. The general assumption for modeling fluid flow through a resting rock matrix (or through a matrix with a much slower velocity than the fluid) is thermal equilibrium, that is, temperature may vary spatially and temporally, but is equal in both phases. For hydrothermal circulation in the crust, the solid rock matrix is usually assumed to be undeformable and at rest. In contrast, for melt segregation in the Earth's crust or mantle the rock matrix is commonly described by a deformable viscous or viscoplastic material. In this case, one needs to address the fluid dynamics of a medium consisting of two immiscible fluids with different flow velocity fields. Such a problem is commonly referred to as two-phase flow. While the general theory for hydrothermal fluid flow in a resting matrix is well developed, the theory of two-phase flow is still in the focus of recent research. McKenzie (1984) and Spiegelman (1993) developed the theory of two-phase flow in the limit of high fluid-matrix viscosity contrasts between melt and matrix and low melt fractions. Bercovici et al. (2001) extended the theory to the case of large fluid fractions. Rudge et al. (2011) included chemical disequilibrium and melting. All these formulations assumed local thermal equilibrium based on the assumption that the thermal diffusion time on the grain scale, or more general, on the scale of a characteristic fluid phase length, is much shorter than the timescale of macroscopic flow and heat transport. However, there might exist cases in which the assumption of local thermal equilibrium is violated. Such cases might be rapid channeling of melt (e.g. Stevenson 1989; Golabek et al. 2008; Katz 2008), and the formation of melt-filled dykes in the crust (Rubin 1998). Also, water-filled fractures in the upper crust, important for hydrothermal fluid circulation, might be associated with diffusive length scales which violate the thermal equilibrium criteria. While thermal equilibrium on the scale of connected pores might be a valid assumption, it is certainty violated on the scale of fluid-filled dykes or fractures for non-isothermal systems. In such cases, hot fluids might intrude into the cooler rock matrix over large distances and exchange heat over some interval in time. Particularity, hot melt will intrude into a rock matrix even if the temperature of the rock matrix is below the melting temperature of the fluid material. The intrusion will continue as long as the temperature of the fluid is above its solidus. In the case of hydrothermal equilibrium fluid flow, the characteristic length of fluid transport in the medium can be related to the grain (or pore) size (Fig. 1a), in the case of non-hydrothermal equilibrium the characteristic length would be a characteristic host rock width ds (Fig. 1b), which is the averaged distance between adjacent solid–fluid interfaces through the solid phase. While the porosity φ might be the same in both cases, the conditions for heat transport differ. Figure 1. View largeDownload slide Sketch to illustrate the validity of the assumption of thermal equilibrium in porous flow: (a) fluid flow through pores and (b) fluid flow through dykes or fractures. Figure 1. View largeDownload slide Sketch to illustrate the validity of the assumption of thermal equilibrium in porous flow: (a) fluid flow through pores and (b) fluid flow through dykes or fractures. Amiri & Vafai (1994), Minkowycz et al. (1999), Nield & Bejan (2006), de Lemos (2016) and others presented approximate formulations for thermal non-equilibrium fluid flow considering only the actual, instantaneous temperature difference between fluid and solid phases for heat conduction between the phases and invoking appropriate heat transfer coefficients. In this paper, we follow a more fundamental physical approach and consider not only the actual temperature difference between the two phases, but include the time-dependent effects due to differential movement of the fluid and the solid (e.g. by melt migration) and related earlier temperature changes and heat exchange between the phases. In other words, we will explicitly account for the full time-dependent history of temperature differences between fluid and solid. This history is microscopically hidden within time-dependent thermal boundary layers with microscopically varying temperature gradients normal to the solid–fluid interfaces. Such gradients will be captured by our time-integration approach. Our theory is applicable to both, fluid flow and melt segregation, in a resting undeformable matrix, and to fluid flow or melt segregation in a deformable matrix. 2 OUTLINE OF THE THEORY In our general formulation of two-phase flow with local thermal non-equilibrium, we do not make a physical distinction between the two moving phases, thus the formulation will be phase symmetric. However, for semantic reasons we will call one phase ‘fluid’ with the subscript f, and the other phase ‘solid’ with the subscript s. The phases may differ in their physical properties, particularly in their densities, viscosities and thermal properties. The volume fraction occupied by the fluid will be addressed by the term ‘porosity’, though it is not necessarily a porosity in its strict sense. We distinguish between a microscopic and a macroscopic scale. On the microscopic scale, the two phases, fluid and solid, have distinct geometries and are separated by interfaces, and the temperature in each phase may vary spatially. However, in our approach, we do not treat this interface explicitly, but, on a macroscopic scale, we regard the two phases as continua which exchange heat (see Fig. 2). Thus, the temperature in each phase may vary spatially and may be different in each phase. On the macroscopic scale, the amount of heat exchange depends on the previous history of the temperature difference between solid and fluid as well as on variables such as the actual local temperature difference, thermal conductivities of the phases and a geometric factor characterizing the ratio of the interface surface to the volume. Figure 2. View largeDownload slide Sketch to illustrate the difference between (a) microscopic and (b) macroscopic scales: (a) distinct fluid-filled structures are resolved and the temperature in each phase might be different and (b) a larger region is considered, where fluid is only present in the central square box (grey pattern) with volumetric fraction φ. Both phases are considered as continua which locally can exchange heat and where the temperatures in the phases can be different. Figure 2. View largeDownload slide Sketch to illustrate the difference between (a) microscopic and (b) macroscopic scales: (a) distinct fluid-filled structures are resolved and the temperature in each phase might be different and (b) a larger region is considered, where fluid is only present in the central square box (grey pattern) with volumetric fraction φ. Both phases are considered as continua which locally can exchange heat and where the temperatures in the phases can be different. More specifically, our approach makes use of a time-dependent reference heat exchange function between the solid and fluid ($$q_s^*( t )$$, see below), which has to be determined on the microscopic scale prior to solving the full problem. When solving the full problem on the macroscopic scale by separating the heat equation into one for the solid and one for the fluid, the heat exchange between the two phases will take into account the actual temperature difference between the phases, as well as all previous temperature differences. These have to be convoluted with the reference heat exchange function to account for thermal relaxation. Thus, our approach requires storing the actual and all previous temperatures of both phases for the complete run. Our approach contains some simplifications compared to solving the full microscopic advection-diffusion problem for porous flow: the detailed 3-D microscopic solid and fluid geometries are replaced by equivalent geometries such as layers, spheres or infinite half-spaces. Thus, the detailed geometry of a porous network cannot be captured. No local heat advection in the fluid by differential flow on the pore scale, only the net Darcy flow is accounted for. We assume that the two phases obey the equations of mass and energy conservation (see e.g. Schmeling 2000). In the present approach, we neglect the melting and solidification process, however, this process might be introduced later by invoking the Stefan problem. The velocities of the solid and the fluid phases can be the results of solving the momentum equations for these phases, for example, a Stokes equation for the deformable porous solid and a Darcy equation for the fluid, but the way they are obtained is not of concern for this study. 2.1 Basic equations To understand how heat is exchanged between the fluid and the solid phases, we will first examine the process on the microscopic scale (see Fig. 3a). We assume that at a certain time the locally averaged temperature of the fluid phase in a medium of porosity φ (i.e. the volumetric fraction of the fluid) is given by Tf, that of the solid is Ts, and that of the interface between fluid and solid is Tb. Local averaging is meant as averaging the temperatures over the respective fluid and solid volumes, and the interface, on a length scale larger than the characteristic pore scale (ds in Figs 1 and 2), but smaller than the macroscopic length scale (L in Fig. 2). Conservation of energy of the matrix phase is given by (e.g. Nield & Bejan 2006; de Lemos 2016)   \begin{equation}\left( {1 - \varphi } \right){\rho _s}{c_{ps}}\left( {\frac{{\partial {T_s}}}{{\partial t}} + {{\vec{v}}_s} \cdot \vec{\nabla }{T_s}} \right) = \vec{\nabla }{\rm{\ }} \cdot \left( {\left( {1 - \varphi } \right){\lambda _s}\vec{\nabla }{T_s}} \right) + {Q_b},\end{equation} (1) and that of the fluid phase is   \begin{equation}\varphi {\rho _f}{c_{pf}} \left( {\frac{{\partial {T_f}}}{{\partial t}} + {{\vec{v}}_f} \cdot \vec{\nabla }{T_f}} \right) = \vec{\nabla } \cdot \left( {\varphi {\lambda _f}\vec{\nabla }{T_f}} \right) - {Q_b},\end{equation} (2) Figure 3. View largeDownload slide (a) Microscopic sketch of the thermal non-equilibrium two-phase flow. The temperatures Ts, Tf, Tband the velocities are variable on this microscopic scale, but constant on the macroscopic scale. (b) Temporal evolution of the local temperature field (shown in the grey ellipsis in Fig. 3a)as due to a change in the macroscopic temperature difference between solid and fluid, as a function of the distance x from the interface. Thick black line for time t = 0 and thinner black lines for successive times. Figure 3. View largeDownload slide (a) Microscopic sketch of the thermal non-equilibrium two-phase flow. The temperatures Ts, Tf, Tband the velocities are variable on this microscopic scale, but constant on the macroscopic scale. (b) Temporal evolution of the local temperature field (shown in the grey ellipsis in Fig. 3a)as due to a change in the macroscopic temperature difference between solid and fluid, as a function of the distance x from the interface. Thick black line for time t = 0 and thinner black lines for successive times. where Qb is the interphase heat exchange (IPHE) rate (i.e. the thermal coupling term of the fluid and solid phase across the interfaces) with the dimension Jm−3 s−1; all other variables are explained in Table 1. In our formulation of eqs (1) and (2), we neglect heat generation due to internal heating, latent heat and dissipation (but see Rudge et al. (2011) for a more complete formulation for the thermal equilibrium case). Note, that in non-equilibrium two-phase flow, we need two separate equations for heat transport. Table 1. Symbols, their definitions and physical units used in this study. Symbol  Definition  Units  as  Geometry factor depending on the dimension of the problem  –  cps  Specific heat capacity of the solid phase  J kg−1 K−1  cpf  Specific heat capacity of the fluid phase  J kg−1 K−1  d  Grain size  M  df  Characteristic length of the fluid phase  M  ds  Characteristic length of the solid phase  M  L  Macroscopic length scale  M  l0  Wavelength for fracture/dyke spacing  M  Pe  Peclet number  –  PeD  Darcy flow related Peclet number  –  PeL  Macroscopic Peclet number  –  Qb  Interphase heat exchange rate (IPHE)  J s−1 m−3  Δq  Heat flow through the fluid–solid interface  W m−2  qs  Specific heat flow from the fluid to the solid phase, normal to the interface  W m−2  $$q_s^*$$  Specific heat flow response function to a temperature step function  W m−2 K−1  qf  Specific heat flow from the solid to the fluid phase, normal to the interface  W m−2  $${q_{{s_{{\rm{inst}}}}}}$$  Instantaneous specific heat flow from the solid to the fluid phase  W m−2  S  Geometry term, describing the interface area per volume of the solid  m−1  Sb  Surface area of all fluid–solid interfaces in a volume  m2  T  Temperature  °C  Tb, Tb0  Temperature and initial temperature of the fluid–solid interface  °C  Tf, Tf0  Temperature and initial temperature of the fluid phase  °C  Ts, Ts0  Temperature and initial temperature of the solid phase  °C  t, t0  Time and initial time  S  tov  ‘Overtake time’ when advective propagation prevails conductive  S  Vg  Volume containing fluid-filled fractures/dykes (entire volume)  m3  vD  Darcy velocity  m s−1  Vs  Volume of the solid phase  m3  $${\vec{v}_f}$$  Velocity of the fluid phase  m s−1  $${\vec{v}_s}$$  Velocity of the solid phase  m s−1  vdiff  Velocity of conductively moving temperature front  m s−1  W  Heat content per volume  J m−3  Ws, Wf  Heat content per volume of the solid phase and of the fluid phase  J m−3  $$\vec{x}$$  Position vector  M  x  Space coordinate normal to the fractures/dykes  M  xinter  x position of the fluid–solid interface  M  z  Space coordinate in direction of the fractures/dykes  M  δM  thickness of the thermal boundary layer  M  δT  Small temperature increment  K  ΔT0  Initial temperature difference between the fluid and the solid phase  K  κs, κf  Thermal diffusivity of the solid phase and of the fluid phase  m2 s−1  λs, λf  Thermal conductivity of the solid phase and the fluid phase  W m−1 K−1  $$\bar{\lambda }$$  Mean thermal conductivity  W m−1 K−1  ρs, ρf  Density of the solid phase and of the fluid phase  kg m−3  φ  Porosity  –  $${\vec{\psi }_s}$$  Particle path in the solid phase  M  $${\vec{\psi }_{s0}}$$  Initial starting point at the flow path  M  $${\vec{\psi }_f}$$  Particle path in the fluid phase  M  Symbol  Definition  Units  as  Geometry factor depending on the dimension of the problem  –  cps  Specific heat capacity of the solid phase  J kg−1 K−1  cpf  Specific heat capacity of the fluid phase  J kg−1 K−1  d  Grain size  M  df  Characteristic length of the fluid phase  M  ds  Characteristic length of the solid phase  M  L  Macroscopic length scale  M  l0  Wavelength for fracture/dyke spacing  M  Pe  Peclet number  –  PeD  Darcy flow related Peclet number  –  PeL  Macroscopic Peclet number  –  Qb  Interphase heat exchange rate (IPHE)  J s−1 m−3  Δq  Heat flow through the fluid–solid interface  W m−2  qs  Specific heat flow from the fluid to the solid phase, normal to the interface  W m−2  $$q_s^*$$  Specific heat flow response function to a temperature step function  W m−2 K−1  qf  Specific heat flow from the solid to the fluid phase, normal to the interface  W m−2  $${q_{{s_{{\rm{inst}}}}}}$$  Instantaneous specific heat flow from the solid to the fluid phase  W m−2  S  Geometry term, describing the interface area per volume of the solid  m−1  Sb  Surface area of all fluid–solid interfaces in a volume  m2  T  Temperature  °C  Tb, Tb0  Temperature and initial temperature of the fluid–solid interface  °C  Tf, Tf0  Temperature and initial temperature of the fluid phase  °C  Ts, Ts0  Temperature and initial temperature of the solid phase  °C  t, t0  Time and initial time  S  tov  ‘Overtake time’ when advective propagation prevails conductive  S  Vg  Volume containing fluid-filled fractures/dykes (entire volume)  m3  vD  Darcy velocity  m s−1  Vs  Volume of the solid phase  m3  $${\vec{v}_f}$$  Velocity of the fluid phase  m s−1  $${\vec{v}_s}$$  Velocity of the solid phase  m s−1  vdiff  Velocity of conductively moving temperature front  m s−1  W  Heat content per volume  J m−3  Ws, Wf  Heat content per volume of the solid phase and of the fluid phase  J m−3  $$\vec{x}$$  Position vector  M  x  Space coordinate normal to the fractures/dykes  M  xinter  x position of the fluid–solid interface  M  z  Space coordinate in direction of the fractures/dykes  M  δM  thickness of the thermal boundary layer  M  δT  Small temperature increment  K  ΔT0  Initial temperature difference between the fluid and the solid phase  K  κs, κf  Thermal diffusivity of the solid phase and of the fluid phase  m2 s−1  λs, λf  Thermal conductivity of the solid phase and the fluid phase  W m−1 K−1  $$\bar{\lambda }$$  Mean thermal conductivity  W m−1 K−1  ρs, ρf  Density of the solid phase and of the fluid phase  kg m−3  φ  Porosity  –  $${\vec{\psi }_s}$$  Particle path in the solid phase  M  $${\vec{\psi }_{s0}}$$  Initial starting point at the flow path  M  $${\vec{\psi }_f}$$  Particle path in the fluid phase  M  View Large The main problem is the derivation of the IPHE term. To determine an expression for the IPHE, we start on the microscopic scale by considering the time-dependent evolution of the temperature field on either side of a fluid–solid interface with the (hot) fluid phase moving with a differential velocity compared to the solid phase (see Fig. 3a). The temperature difference between solid and fluid may also be caused by different conductive heat fluxes through the two phases due to different thermal parameters. The principal assumption is that at a certain time t0 an initial temperature difference ΔT exists between the solid and fluid phases (e.g. due to hot melt intrusion). This temperature difference initiates a temporally and spatially varying temperature field in either phase on the microscopic scale. For successive times (t > t0), the initial temperature difference will equilibrate (see Fig. 3b). Thus, any temperature difference between a solid and a fluid phase generally represents a transient stage. If we assume a positive temperature change in a hot fluid or melt, it acts as a heat source in the solid (eq. 1) and a sink in the fluid (eq. 2, i.e. the hot fluid is losing heat to the colder solid). It depends on the local heat flux across the interface and the specific interfacial area S of dimension m−1 describing the interface area per volume   \begin{equation}{Q_b} \left( {\vec{x},t} \right) = \ S \cdot {q_s}\left( {\vec{x},t} \right).\end{equation} (3) qs is the heat flow normal to the interface from the fluid phase into the solid phase (see Fig. 3a). Since we only consider the normal component of the heat flow it is treated in the following as a scalar variable. Qb could be determined in a similar way by using the interface normal of the heat flow from the solid phase to the fluid phase qf. We will first derive an expression for $${q_s}( {\vec{x},t} )$$ and further down discuss the specific interfacial area term S. 2.2 Interface heat flux In the following, we treat the heat flux resulting from a sudden temperature change as described by the ‘heat flow response function to a temperature step function’, denoted by an asterisk $$q_s^*( t )$$  \begin{equation}{q_s} \left( t \right) = \ \Delta {T_0} \cdot q_s^*\left( t \right).\end{equation} (4) The expression for $$q_s^*( t )$$ depends on the application, i.e. the choice of the phase geometries. A possible choice could be the assumption of a cooling half-space ($$q_s^*( t ) = {\lambda _s}\frac{1}{{\sqrt {\pi {\kappa _s}t} }}$$, e.g. Turcotte & Schubert 2002). Another possibility is the assumption of a series of cooling (vertical or horizontal) sheets (for which we will derive the expression further down). As a result of an initial temperature difference  ΔT0 at all interfaces in a volume, a subsequent heat flux from the fluid into the solid,$${\rm{\ \Delta }}{T_0}q_{s0}^*( {\tilde{t}} )$$, is generated, where $$\tilde{t}$$ is the subsequent time starting with t0, that is, $$\tilde{t} = \ t - {t_0}$$. For subsequent times  t > t0, we determine the heat flux into the solid by considering any temporal changes δT(τ) of the macroscopic fluid–solid temperature difference along a solid volume flow path (see Fig. 4) taking place at any time in the history of the solid particle τ < t. Each of these incremental temperature changes cause an additional time-dependent heat flow into the solid, which, at the time t, has decayed to $$d{q_{s,\tau }}( t ) = \ \delta T( \tau )q_s^*( {t - \tau } )$$. The full interface heat flux at time t can thus be determined by summing up $${\rm{\Delta }}{T_0}q_{s0}^*( {\tilde{t}} )$$ and an integral over all dqs, τ(t) having emerged at all times τ < t. The problem is that this integration has to take into account the history of solid and fluid flow paths, and that the temperature changes δT(τ) need some further attention. Figure 4. View largeDownload slide Solid volume flow path $${\vec{\psi }_s}( \tau )$$and fluid volume flow paths $${\vec{\psi }_{fi}}( \tau )$$, where the fluid flow paths cross the solid flow path at different times. The fluid flow path $${\vec{\psi }_{fn}}( \tau )$$merges the solid path at $$( {\vec{x},t} )$$. At previous time τi , the solid particle has been exposed to different fluid particles, denoted by the fluid flow paths. At those times, the solid particle has experienced incremental temperature changes $${ {\delta T} |_{{{\vec{\psi }}_s}( \tau )}}$$ whose subsequent response needs to be determined along the solid flow path all the way to $$( {\vec{x},t} )$$. Figure 4. View largeDownload slide Solid volume flow path $${\vec{\psi }_s}( \tau )$$and fluid volume flow paths $${\vec{\psi }_{fi}}( \tau )$$, where the fluid flow paths cross the solid flow path at different times. The fluid flow path $${\vec{\psi }_{fn}}( \tau )$$merges the solid path at $$( {\vec{x},t} )$$. At previous time τi , the solid particle has been exposed to different fluid particles, denoted by the fluid flow paths. At those times, the solid particle has experienced incremental temperature changes $${ {\delta T} |_{{{\vec{\psi }}_s}( \tau )}}$$ whose subsequent response needs to be determined along the solid flow path all the way to $$( {\vec{x},t} )$$. We consider that the fluid and possibly also the solid move on independent paths (in a fixed Eulerian coordinate system). The heat fluxes and temperatures discussed above, however, are related to the heat contents of fluid and solid volumes, which will be carried away by either flow. The position of a moving solid volume particle can be described by a flow path (time-dependent position vector) $${\vec{\psi }_s}( t )$$ (see Fig. 4), starting at an initial position $${\vec{\psi }_{s0}}$$. We have to consider all incremental fluid–solid temperature changes (δT = δTf − δTs) as incremental step functions along such a solid particle volume flow path, which are not associated with the instantaneous interphase heat exchange (IPHE), i.e. the heat exchange along the solid flow line at a particular time (denoted further down by τ). In other words, δTs and δTf may be due to some ‘external’ macroscopic heat transport effects such as advection, conduction, heat generation and phase changes. In particular, we need to consider processes which act differently in the two phases, for example, advection with different $${\vec{v}_s}$$ and $${\vec{v}_f}$$, conduction with different thermal conductivities. For the incremental temperature change δT along the path at times τ ≤ t, we can write   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} - \left( {{{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}}} \right)} \right)\ \delta \tau ,\end{equation} (5) where $${ {\frac{{d{T_f}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ are the entire temperature changes in the solid and the fluid phases along the path $${\vec{\psi }_s}( \tau )$$ (i.e. following the solid particle), respectively, and $${ {\frac{{d{T_f}}}{{dt}}} |_{{Q_b},{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{Q_b},{{\vec{\psi }}_s}( \tau )}}$$ are the changes due to the instantaneous interphase heat exchange. While $${ {\frac{{d{T_f}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ and $${ {\frac{{d{T_s}}}{{dt}}} |_{{{\vec{\psi }}_s}( \tau )}}$$ can simply be evaluated along the flow path, an expression for $$( { - {{ {\frac{{d{T_f}}}{{dt}}} |}_{{Q_b},{{\vec{\psi }}_s}( \tau )}} + {{ {\frac{{d{T_s}}}{{dt}}} |}_{{Q_b},{{\vec{\psi }}_s}( \tau )}}} )$$ needs to be derived. The heat contents W per volume (including both phases) at any instant of time is given for the solid and the fluid phases by   \begin{equation}{W_s} = \left( {1 - \varphi } \right)\ {\rho _s}{c_{Ps}}{T_s} \,{\rm{and}}\,{W_f} = \varphi {\rho _f}\ {c_{Pf}}{T_f}.\end{equation} (6) These heat contents change locally with time, and the contribution due to the instantaneous interphase heat exchange rate (denoted by subscript Qb) can written in terms of the IPHE (eq. 3) as   \begin{equation}S \cdot \ {q_s} = {\left. {\frac{{d{W_s}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = \ - \ {\left. {\frac{{d{W_f}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {1 - \varphi } \right)\ {\rho _s}{c_{Ps}}\ {\left. {\frac{{d{T_s}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}} = - \varphi {\rho _f}{c_{Pf}}{\left. {\frac{{d{T_f}}}{{dt}}} \right|_{{Q_b},{{\vec{\psi }}_s}\left( \tau \right)}}.\end{equation} (7) Rearranging for the temperature change rates and inserting them in eq. (5) leads to   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {\frac{{S \cdot {q_s}}}{{\varphi {\rho _f}{c_{Pf}}}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {\frac{{S \cdot {q_s}}}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}}} \right)\delta \tau \,,\end{equation} (8) using again eq. (3), this is   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = \left( {{{\left. {\frac{{d{T_f}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} - {{\left. {\frac{{d{T_s}}}{{dt}}} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}} + {{\left. {{Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right|}_{{{\vec{\psi }}_s}\left( \tau \right)}}} \right)\ \delta \tau .\end{equation} (9) To consider all incremental temperature changes δT along a solid volume path, we sum up the initial heat flow response, $${\rm{\Delta }}{T_0}q_{s0}^*( {t - {t_0}} )$$, and the $$d{q_{s,\tau }}( t ) = \ \delta T( \tau )q_s^*( {t - \tau } )$$ increments integrated over all previous times τ. This integral, $$\mathop \smallint \nolimits_{{t_0}}^t q_s^*( {t - \tau } ){ {\frac{{\delta T}}{{\delta \tau }}} |_{{{\vec{\psi }}_s}( \tau )}}\,d\tau $$, has to account for the temperature changes δT(τ) along the solid volume path $${\vec{\psi }_s}( \tau )$$ from the initial time t0 to the actual time t. Taking   \begin{equation}{\left. {\delta T} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} = {\left. {\frac{{\partial T}}{{\partial t}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot {\left. {\frac{{dt}}{{d{{\vec{\psi }}_s}}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot d{\psi _s} = {\left. {\frac{{\partial T}}{{\partial t}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot {\left. {\frac{1}{{{{\vec{v}}_s}}}} \right|_{{{\vec{\psi }}_s}\left( \tau \right)}} \cdot d{\psi _s},\end{equation} (10) where $${\vec{v}_s}$$ is the velocity of the solid phase, and using eq. (9), the sum of the initial heat flow response and the integrated heat flow increments yields   \begin{equation}{q_s}\!\left( {{{\vec{\psi }}_s}\left( t \right)} \right) = \ \Delta {T_0}\ q_{s0}^*{\left( {t - {t_0}} \right)_{{{\vec{\psi }}_{s0}}}} + \int_{{{{\vec{\psi }}_{s0}}}}^{{{{\vec{\psi }}_s}\left( t \right)}}{{q_s^*\left( {t - \tau \left( {\psi _s^\prime} \right)} \right)}} \cdot \ {\left( {\frac{{d{T_f}}}{{dt}} - \frac{{d{T_s}}}{{dt}} + {Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right)_{\psi _s^{'}}} \cdot {\left. {\frac{1}{{{{\vec{v}}_s}}}} \right|_{\psi _s^{'}}}d\psi _s^{'}.\end{equation} (11)$$\tau ( {{{\vec{\psi }}_s}} )$$ is the time at which the considered solid particle volume passed a previous position $${\vec{\psi }_s}$$. In case of a resting solid and only a moving fluid, eq. (11) degenerates to a convolution integral   \begin{equation}{q_s}\!\left( {\vec{x},t} \right) = \ \Delta {T_0}\ q_{s0}^*\left( {t - {t_0}} \right) + \int_{{{t_0}}}^{t}{{q_s^*\left( {t - \tau } \right)}} \cdot \ {\left( {\frac{{d{T_f}}}{{dt}} - \frac{{d{T_s}}}{{dt}} + {Q_b}\left( {\frac{1}{{\varphi {\rho _f}{c_{Pf}}}} + \frac{1}{{\left( {1 - \varphi } \right){\rho _s}{c_{Ps}}}}} \right)} \right)_\tau }d\tau ,\end{equation} (12) eq. (12) allows determining the interphase heat exchange rate at a given position and time by integrating over all previous temperature changes of the fluid, which has passed that position, as well as integrating the history of solid temperature changes and heat exchange at that position. 2.3 Heat exchange for a medium containing fluid layers When discussing eq. (4), we already mentioned the choice of a cooling half-space solution for $$q_s^*( t )$$. This is a suitable choice if the underlying assumption is related to widely spaced fluid channels, dykes or fractures (or fluid layers in general) with heat exchange to the solid but no heat exchange in between the fractures. For densely spaced fluid layers such an assumption does not hold. Here, we develop an expression for heat flow from the fluid into the solid for a case of densely spaced fluid layers (or vertical dykes) where the layer spacing influences the overall heat exchange. If we recall Fig. 3(b), a case of a swarm of large vertical fluid layers can be constructed as follows: since the layers are large in vertical direction and we regard only the normal component of the heat flow from the fluid-filled layers into the solid, we can treat the problem as 1D. For a swarm of fluid layers of equal spacing ds and width df, the initial or any incremental temperature difference between fluid and solid which is advected to a certain macroscopic position can be regarded as an initial box-shaped temperature function on microscopical scale. The related initial temperature as a function of a space coordinate x normal to the fluid layer is described by the Fourier series of a rectangular half-wave with an initial temperature difference between fluid and solid phases of ΔT0 and fluid layer spacing given by a wavelength l0 (see Fig. 5)   \begin{equation}{T_0} \left( x \right) = \frac{{{d_f} \cdot \Delta {T_0}}}{{{l_0}}}\ + \mathop \sum \limits_{n\ = \ 1}^\infty {\frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \frac{{{d_f}n\pi }}{{{l_0}}}\cos \frac{{2\pi nx}}{{{l_0}}}} .\end{equation} (13) Figure 5. View largeDownload slide Sketch to illustrate the assumption of a periodic temperature perturbation to model narrowly spaced fluid layers such as dykes or fractures filled with a hot fluid. Figure 5. View largeDownload slide Sketch to illustrate the assumption of a periodic temperature perturbation to model narrowly spaced fluid layers such as dykes or fractures filled with a hot fluid. Since the macroscopic porosity can be expressed by $$\varphi \ = \frac{{{d_f}}}{{{l_0}}}$$, eq. (13) is rewritten to   \begin{equation}{T_0} \left( x \right) = \ \Delta {T_0}\varphi + \mathop \sum \limits_{n\ = \ 1}^\infty {\frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \left( {n\pi \varphi } \right)\cos \frac{{2\pi n\varphi x}}{{{d_f}}}.} \end{equation} (14) Assuming the same thermal properties in the solid and fluid, the temperature for successive times t is determined by solving the 1D conductive heat equation $$\frac{{\partial T}}{{\partial t}} = {\kappa _s}\frac{{{\partial ^2}T}}{{\partial {x^2}}}$$ with the initial condition (eq. 14) for each Fourier term separately. The result is Tn (x, t) = T0n (x) · exp ( − bnt) where T0n(x) is the Fourier terms of the initial condition in eq. (14) and the coefficients bn are obtained as $${b_n} = {\kappa _s}{( {\frac{{2\pi n\varphi }}{{{d_f}}}} )^2}$$. Thus, the temperature in a continuum containing initially rectangular half-wave temperature differences is   \begin{equation}T ( {x,t} ) = \Delta {T_0}\varphi + \mathop \sum \limits_{n\ = \ 1}^\infty \frac{{2 \cdot \Delta {T_0}}}{{n\pi }}\sin \left( {n\pi \varphi } \right)\cos \frac{{2\pi n\varphi x}}{{{d_f}}}\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (15) Applying Fourier's law of heat conduction in one dimension $$q\ = \ - {\lambda _s}\frac{{dT}}{{dx}}$$ to eq. (15) gives the heat flow in such a continuum   \begin{equation}q ( {x,t} ) = {\lambda _s}\ \frac{{4\Delta {T_0}\varphi }}{{{d_f}}}\mathop \sum \limits_{n\ = \ 1}^\infty \sin \left( {n\pi \varphi } \right)\sin \frac{{2\pi n\varphi x}}{{{d_f}}}\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (16) Since we are interested in the heat flow qs from the fluid phase into the solid phase at the interface of the two phases, which is defined by $$x\ = {x_{{\rm{inter}}}}\ = \frac{{{d_f}}}{2}$$ (see Fig. 5), we evaluate eq. (16) at x = xinter  \begin{equation}{q_{s\_{\rm{Layer}}}}\ \left( t \right) = {\lambda _s}\ \frac{{4\Delta {T_0}\varphi }}{{{d_f}}}\mathop \sum \limits_{n\ = \ 1}^\infty {\sin ^2}\left( {n\pi \varphi } \right)\exp \left( { - {\kappa _s}t{{\left( {\frac{{2\pi n\varphi }}{{{d_f}}}} \right)}^2}} \right).\end{equation} (17) This is the exchange heat flux of the layered solid-fluid case. Eq. (17) can be written in the form of   \begin{equation}q_s^* \left(t \right) = \frac{1}{{\Delta {T_0}}}\ {q_{s\_{\rm{Layer}}}}\left( t \right)\end{equation} (18) (see eq. 4) to be used as the ‘heat flow response function to a temperature step function’ in eq. (11) and in eq. (12). 2.4 The geometry term The geometry term S in eq. (3) has the dimension m−1 and is defined by the surface area Sb of all fluid–solid interfaces in a volume element to the entire volume Vg  \begin{equation}S\ = \frac{{{S_b}}}{{{V_g}}}.\end{equation} (19) For spherical solid inclusions of diameter ϑ (i.e. radius ϑ/2), the interface surface is given by Sb = πϑ2 and the solid phase volume by $${V_s} = \frac{1}{6}\ \pi {\vartheta ^3}$$. Since Vg = Vs(1 − φ)− 1, S for spherical solid inclusions is S = 6ϑ− 1(1 − φ). One can do this little exercise for other geometries to find the general relation   \begin{equation}S\ = {a_s}\ d_s^{ - 1}\left( {1 - \varphi } \right),\end{equation} (20) where, more generally, ds is a scaling length of the solid phase of the composite material (i.e. the average distance between adjacent solid–fluid interfaces through the solid phase) and as ≈ 6 for 3-D (e.g. spherical solid volumes), as ≈ 4 for 2-D (e.g. cylindrical solid volumes) and as ≈ 2 for 1-D (e.g. solid and fluid layers). Our simplifying assumption for deriving eq. (20) implies complete wetness, that is, solid grains (or inclusions) are totally surrounded by fluid films. For the 2-D case of parallel fluid dykes, which we study here, this assumption is approximately valid. For other melt geometries such as melt pockets or tubes at the edges of solid grains, the wetness is less than 1 and a different formula than eq. (20) would emerge. 2.5 Non-dimensionalization of the equations For using the equations above in a numerical code and studying the principal behaviour of non-equilibrium two-phase flow, it is convenient to use the equations in non-dimensionalized form. We do this by choosing a set of (dimensional) scaling parameters ($$\skew4\tilde{d},\ \skew4\tilde{\kappa },\ \skew4\tilde{\rho },\ {\skew4\tilde{c}_p},\ \skew4\tilde{T}$$) and setting   \begin{equation}x = x{'} \cdot \tilde{d},\qquad t = t{'} \cdot \frac{{{{\tilde{d}}^2}}}{{\tilde{\kappa }}},\qquad \vec{v} = \vec{v}{'} \cdot \frac{{\tilde{\kappa }}}{{\tilde{d}}},\qquad T = T{'} \cdot \tilde{T},\qquad \lambda = \lambda {'} \cdot \tilde{\kappa }\,\tilde{\rho }\,{\tilde{c}_p},\end{equation} (21) where the apostrophe indicates non-dimensional variables. By specifying $$\tilde{d} = {d_s}\ ,\ \ \widetilde {\kappa \ } = {\kappa _s}\ ,\ \ \tilde{\rho } = {\rho _s}\ ,\ \ {\tilde{c}_p} = {c_{ps}}\ ,\ \ \tilde{T} = \ {\rm{\Delta \ }}{T_0} = {T_{f0}}\ - {T_{s0}}$$, we scale the initial temperature difference for the fluid phase and the solid phase to 1 and obtain for the heat flow, the specific interface S, and Qb  \begin{equation}q = {\kappa _s}\,{\rho _s}\,{c_{ps}}\frac{{\Delta {T_0}}}{{{d_s}}}q{'},\quad S = \frac{1}{{{d_s}}}S{'},\quad {Q_b} = {\kappa _s}\,{\rho _s}\,{c_{ps}}\frac{{\Delta {T_0}}}{{{d_s}^2}}{Q{'}_b}.\end{equation} (22) Non-dimensionalizing of eq. (18) and considering that $${d_f} = \frac{\varphi }{{( {1 - \varphi } )}}{d_s}$$ leads to   \begin{equation}q{_s^{*\prime} } \left( {t{'}} \right) = 4\left( {1 - \varphi } \right)\mathop \sum \limits_{n\ = \ 1}^\infty {\sin ^2}\left( {n\pi \varphi } \right)\exp \left( { - t{'}{{\left( {2\pi n\left( {1 - \varphi } \right)} \right)}^2}} \right).\end{equation} (23) We also need non-dimensionalized versions of eq. (1) for the solid phase   \begin{equation}\left( {1 - \varphi } \right)\ \left( {\frac{{\partial {T_s}^\prime }}{{\partial t{'}}} + {{\vec{v}}_s}{'} \cdot \vec{\nabla }{'}{T_s}{'}} \right) = \vec{\nabla }{'} \cdot \left( {\left( {1 - \varphi } \right)\vec{\nabla }{'}{T_s}{'}} \right) + {Q_b}^\prime ,\end{equation} (24) and eq. (2) for the fluid phase   \begin{equation}\varphi \ \left( {\frac{{\partial {T_f}{'}}}{{\partial t{'}}} + \vec{v}_f^{'} \cdot \vec{\nabla }{'}{T_f}{'}} \right) = \vec{\nabla }{'} \cdot \left( {\varphi \vec{\nabla }{'}{T_f}^\prime } \right) - {Q_b}^\prime .\end{equation} (25) The non-dimensionalized form of the geometry term is S' = as (1 − φ), and eq. (12) with the convolution integral is given by   \begin{equation}{q_s}^\prime \ \left( {\vec{x}{'},t{'}} \right) = q{_{s0}^{*\prime} }\left( {t{'} - {{t{'}}_0}} \right) + \mathop \smallint \nolimits_{{t_0}{'}}^{t{'}} q_s^*{'}\left( {t{'} - \tau {'}} \right) \cdot \ {\left( {\frac{{d{T_f}{'}}}{{dt{'}}} - \frac{{d{T_s}{'}}}{{dt{'}}} + {Q_b}{'}\left( {\frac{1}{{\varphi {\rho _f}{'}{c_{Pf}}{'}}} + \frac{1}{{\left( {1 - \varphi } \right)}}} \right)} \right)_{\tau {'}}}d\tau {'}.\end{equation} (26) 3. APPLICATION OF THE THEORY First, we determined the maximum Fourier coefficient nmax needed in the sum in eq. (23) to get a sufficient approximation for the response function for the case of parallel fluid layers (cf. Fig. 5). We tested values between 10 and 500 and found that >100 is necessary for reliable results. In this study, we used nmax = 500. Fig. 6 shows the non-dimensional heat flow response function $$q_s^{*{'}}$$ as a function of time and porosity in a log–log plot. Typically, three or two stages can be identified for all curves in Fig. 6, characterized by different slopes in the $$\log q_s^{*{'}} - \log t{'}$$ plot. The first stage, independent of the porosity, has a slope of $${\sqrt {t{'}} ^{ - 1}}$$ and represents the solution of two half-spaces of different initial temperature brought into sudden contact. Heat exchange is not yet influenced by the neighbouring fluid layer and the central regions in the solid and fluid layers are still at their initial temperatures. The second stage has a slope of about (t')− 3/2 (this can be verified by choosing porosities closer to 0 or 1) and reflects the stage in which the thermal boundary layers in the thinner layer (fluid or solid, depending on φ) have merged in the centre of the thinner layer, while the centre of the thicker layer is still not affected by the cooling/heating. In contrast to the first stage, the temperature at the boundary between the fluid and solid drops (or rises) with time, which leads to a stronger decay than proportional to $${\sqrt {t{'}} ^{ - 1}}$$. The last stage has a slope according to e− t'and describes the situation when heat flux rapidly ceases. For equally thick solid and fluid layers (φ = 0.5), the second stage is skipped, because the thermal boundary layers merge in the centres of both layers simultaneously. The differences of the curves in respect to porosity is caused by the difference in total heat contents of the fluid layers, since their volume is a function of the porosity while the heat exchange boundary is not in this 1-D case. Figure 6. View largeDownload slide The non-dimensional fluid–solid interface ‘heat flow response function to a temperature step function’ $$q_s^{*{'}}$$ (eq. 23) as a function of the non-dimensional time t'and the volumetric fraction of the fluid φ (‘porosity’). The sum in eq. (23) was truncated at used nmax = 500. The three different stages which are discussed in the text are indicated. Figure 6. View largeDownload slide The non-dimensional fluid–solid interface ‘heat flow response function to a temperature step function’ $$q_s^{*{'}}$$ (eq. 23) as a function of the non-dimensional time t'and the volumetric fraction of the fluid φ (‘porosity’). The sum in eq. (23) was truncated at used nmax = 500. The three different stages which are discussed in the text are indicated. To apply our theory for non-equilibrium two-phase flow, we studied two simple 1-D test cases. For all cases, we used $$q_s^{*{'}}( {t{'}} )$$ as given in eq. (23) and shown in Fig. 6, and solved eq. (26) for the heat flow from the fluid phase to the solid phase. With the result for $$q_s^{'}( {\vec{x},t} )$$, we determined the IPHE term $$Q_b^{'}( {\vec{x},t} )$$ (using S' and $$q_s^{'}( {\vec{x},t} )$$ in eq. 3) and setting as = 2. With this $$Q_b^{'}( {\vec{x},t} )$$, we solved the coupled heat transport equations (eqs 24 and 25) for the fluid and the solid phases, respectively. We studied the following test cases and assumed for the sake of simplicity the same thermal conductivity for both phases, λf = λs: Case 1: solid and fluid phases in rest. This is accomplished by setting $${\vec{v}_f} = \ 0$$ in eq. (25). Case 2: solid phase in rest, fluid phase moving. Case 2a evaluates the full convolution integral $$\mathop \smallint \nolimits_{{t_0}}^t\, (\,\,\,)\,d\tau $$ (eq. 12), case 2b traces the integral only for a fraction of the diffusion time, ds2/κs. For all cases, we compare our result with a formulation (Amiri & Vafai 1994; Minkowycz et al. 1999) in which the heat transfer between fluid and solid is described by the instantaneous temperature difference Tf–Ts and an interstitial heat transfer coefficient. Their formulation can be rewritten as   \begin{equation}{q_{{s_{{\rm{inst}}}}}}\! \left( {\vec{x},t} \right) = \frac{{\bar{\lambda }}}{{{\delta _M}}}\ \left( {{T_f}\left( {\vec{x},t} \right) - {T_s}\left( {\vec{x},t} \right)} \right),\end{equation} (27) where $$\bar{\lambda }$$ is the mean thermal conductivity of the solid and the fluid phases, in general an arithmetic ($$\bar{\lambda } = \ \varphi {\lambda _f} + ( {1 - \varphi } ){\lambda _s}$$) or geometric ($$\bar{\lambda } = \lambda _f^\varphi \ \cdot \lambda _s^{( {1 - \varphi } )}$$) average may be used, and δM is the microscopic thickness of the thermal boundary layer between the fluid and the solid phases. As indicated in eq. (27), the interstitial heat transfer coefficient can be understood as a mean thermal conductivity divided by an effective thermal boundary layer thickness δM. In the formulations by Minkowycz et al. (1999) and Amiri & Vafai (1994), the interstitial heat transfer coefficient is formulated in terms of constant model-dependent non-dimensional numbers such as the Sparrow number, Reynolds number and Prandtl number (for definitions, see their papers). From our formulation, however, it is clear that the effective thermal boundary layer is time dependent. It is thus interesting to test, to which extent constant values of δM may approximate the thermal non-equilibrium behaviour in porous flow. We use the subscript ‘inst’ for the resulting heat flow, since this formula gives the instant response to a temperature difference between the two phases, but does not consider the thermal history of the interphase heat exchange. The equations were programmed in MATLAB (MATLAB, 2012), using second-order finite differences (FDs) for the numerical equivalent of eqs (24) and (25), with upwind for the advection term in eq. (25), various integration rules for eq. (26) and a maximum Fourier coefficient nmax = 500 in eq. (23). Integration of eq. (26) needs special attention because the integrand has a singularity at the upper boundary due to the effective $${\sqrt {t'} ^{ - 1}}$$ behaviour of q*. In principle, for a finite number nmax of Fourier coefficients q* remains finite and the trapezoidal rule may be used to evaluate eq. (26), though with very small time steps (cf.Section 3.1). Alternatively, for case 2, we integrate eq. (26) using the Simpson rule up to the time t' − dt', where dt' is the time step (for an odd number of integration intervals we apply the Trapez rule for the very first interval, because the Simpson rule requires an even number of intervals). Rather than solving the remaining improper integral $$\mathop \smallint \nolimits_{t{'} - dt{'}}^{t{'}} A( {\tau {'}} )q_s^*( {t{'} - \tau {'}} )d\tau {'}$$ numerically we assume $$q_s^*( \tau )\ = \ a/\sqrt \tau $$ and utilize the analytic solution (using the mean value theorem) $$\mathop \smallint \nolimits_0^{d{t^{'}}} A( \tau )q_s^*( \tau )d\tau \ = \ 2A( {{\tau _{{\rm{mean}}}}} )q_s^*( {dt{'}} )dt{'}$$, that is, the argument A has to be taken at an intermediate time 0 ≤ τ'mean ≤ dt' and q*(dt') is taken at the first time increment dt'. In our calculations, we use τ'mean = 0.3dt' but other choices gave similar results. The numerical values of the non-dimensionalized variables are given in Table 2. Table 2. Numerical values of the non-dimensional variables. Symbol  Definition  Value  as  Geometry factor for the solid phase  2  $$d_s^{'}$$  Characteristic length of the solid phase  1  φ  Porosity  0.1–0.9  h'  Width of the model (depending on diffusion time)  0.125–10  $$T_{f0}^{'}$$  Initial temperature of the fluid phase  1  $$T_{s0}^{'}$$  Initial temperature of the solid phase  0  $${\rm{\Delta }}T_0^{'}$$  Initial Temperature difference between the fluid and the solid phases  1  δM'  Effective thermal boundary layer thickness  0.05–3  $$t_{{\rm{max}}}^{'}$$  Maximum time  0.5–3  dt'  Time step  0.1–0.00001  $$v_f^{'}$$  Velocity of the fluid (i.e. Peclet number)  0–300  dx'  Grid size of the (macroscopic) model  0.01–0.02  Symbol  Definition  Value  as  Geometry factor for the solid phase  2  $$d_s^{'}$$  Characteristic length of the solid phase  1  φ  Porosity  0.1–0.9  h'  Width of the model (depending on diffusion time)  0.125–10  $$T_{f0}^{'}$$  Initial temperature of the fluid phase  1  $$T_{s0}^{'}$$  Initial temperature of the solid phase  0  $${\rm{\Delta }}T_0^{'}$$  Initial Temperature difference between the fluid and the solid phases  1  δM'  Effective thermal boundary layer thickness  0.05–3  $$t_{{\rm{max}}}^{'}$$  Maximum time  0.5–3  dt'  Time step  0.1–0.00001  $$v_f^{'}$$  Velocity of the fluid (i.e. Peclet number)  0–300  dx'  Grid size of the (macroscopic) model  0.01–0.02  View Large 3.1 Case 1: solid and fluid phases in rest The first test case is an initial value problem with $$T_f^{'}( {t{'} = \ 0} ) = \ 1$$ and $$T_s^{'}( {t{'} = \ 0} ) = \ 0$$. Both the fluid and the solid are assumed at rest. On the macroscopic scale, no temperature variations are assumed, thus heat eqs (24) and (25) degenerate to   \begin{equation}\left( {1 - \varphi } \right)\ \frac{{\partial {T_s}{'}}}{{\partial t{'}}} = \ + {Q_b}^\prime \end{equation} (28) and   \begin{equation}\varphi \ \frac{{\partial {T_f}{'}}}{{\partial t{'}}} = \ - {Q_b}^\prime .\end{equation} (29) Before solving these equations with our method, we derive the analytical solution by switching to the microscopic scale. On the microscopic scale, the fluid is assumed to occupy infinite planar fluid layers, thus we have the rectangular temperature distribution discussed previously (see eq. 15) and Fig. 5). On the microscopic scale, the temperature in the fluid and the solid phases, respectively, can directly be evaluated by solving the non-dimensional form of eq. (15). To obtain the macroscopic fluid and solid temperatures, the obtained temperatures are averaged by spatial integration of eq. (15) over the fluid and solid volumes, respectively. The resulting temperatures represent the spatially constant fluid and solid temperatures on the macroscopic scale and are functions of time and of porosity (which is in our case the volumetric fraction of the fluid). The result is depicted in Fig. 7 and will be called in the following ‘the analytical solution’. Figure 7. View largeDownload slide Non-dimensional temperature versus non-dimensional time for the fluid phase (red) and the solid phase (blue) as a function of time for different porosities. The porosities φ valid for each curve are depicted at the end of the curves, outside of the box. The initial temperature disturbance is a rectangular half-wave. Figure 7. View largeDownload slide Non-dimensional temperature versus non-dimensional time for the fluid phase (red) and the solid phase (blue) as a function of time for different porosities. The porosities φ valid for each curve are depicted at the end of the curves, outside of the box. The initial temperature disturbance is a rectangular half-wave. For thin fluid layers (φ ∼ ≤ 0.3), the fluid is rapidly cooled, when coming into contact with the solid material, which in turn heats up until thermal equilibrium is reached at times t' ≤ 0.3. On the contrary, for wide fluid layers (φ ∼ ≥ 0.7), the temperature decrease in the fluid is much slower and equilibrium is reached after a much longer time. To test our approach, we calculate the IPHE term with the convolution integral (eq. 26), and solve the heat eqs (28) and (29) for porosities φ equal to 0.1, 0.5 and 0.9 and compare the results to the analytical solution and the approximate solutions based on eq. (27). The results are shown in Fig. 8. For each porosity case, we used three different non-dimensional time steps dt' for the integration of the convolution to show the convergence of our approach. The inset in Fig. 8(c) shows the convergence of the errors of the late equilibrated stages in a log–log diagram. For φ = 0.1 (Fig. 8a), the dykes are very thin and cool rapidly. On the microscopic scale, this behaviour is associated with a very steep temperature gradient in the boundary layer between the fluid and the solid phases. On the macroscopic scale, rapid cooling demands a sufficiently high temporal resolution, here small time steps of about 10−5 are needed, the error order is about 1/2 (Fig. 8a and inset in c). The approximate solution with instantaneous heat exchange (eq. 27) using δM = 0.35 underestimates the thermal equilibration in the fluid and the solid phases at early times (t' < 0.04), but gives good results for later times (Fig. 8d). For φ = 0.5 (Fig. 8b and inset in c), a time step of dt' = 10−4 gives a very good fit to the analytical solution with an error order also of about 1/2; the approximate solution still underestimates the thermal equilibration in the very early cooling stage, but is also very close to the analytical solution. Figure 8. View largeDownload slide Non-dimensional temperature for the solid phase (blue lines) and the fluid phase (red lines) versus non-dimensional time: (a)–(c) our approach: based on an initial rectangular half-wave temperature perturbation and (d)–(f) approximate solution: based on eq. (27), for porosities of (a) and (d) 0.1, (b) and (e) 0.5 and (c) and (f) 0.9 in comparison to the analytical solution (green lines, see Fig. 7). For (a)–(c) three different non-dimensional time steps dt΄ are used for each model to show the convergence of the approach; for (d)–(f) four different δ'M values are used (and the smallest time steps dt΄ according to the results in (a)-(c). Note the difference in t΄ scale. Figure 8. View largeDownload slide Non-dimensional temperature for the solid phase (blue lines) and the fluid phase (red lines) versus non-dimensional time: (a)–(c) our approach: based on an initial rectangular half-wave temperature perturbation and (d)–(f) approximate solution: based on eq. (27), for porosities of (a) and (d) 0.1, (b) and (e) 0.5 and (c) and (f) 0.9 in comparison to the analytical solution (green lines, see Fig. 7). For (a)–(c) three different non-dimensional time steps dt΄ are used for each model to show the convergence of the approach; for (d)–(f) four different δ'M values are used (and the smallest time steps dt΄ according to the results in (a)-(c). Note the difference in t΄ scale. For φ = 0.9 (Fig. 8c), the wide hot dykes heat up the thin solid phase in between the dykes. For this case, the convergence is poor (inset in Fig. 8c). Comparing Fig. 8(c) to (a), equilibration appears to be slow. However, since $$t{'} = \ t \cdot \frac{{\tilde{\kappa }}}{{d_s^2}}$$, t'scales with the quadratic inverse of the characteristic dykes distance ds which is different for the various φ − cases. In fact, the equilibration timescales with the diffusion time associated with the thinner layer, which, in non-dimensional form, is equal 1 in the case of φ = 0.9, but which is equal 0.1 in case of φ = 0.1. Using a time step of dt' = 10−3, we are very close to the analytical solution. Coarser time steps lead to an underestimation of the heat flow and thus, also underestimate the temperature in the solid for early times and lead to an overshot of the temperature at later times. The approximate solution (eq. 27) (Amiri & Vafai 1994; Minkowycz et al. 1999) gives good results for late times, but does not catch the thermal equilibration at early times, neither with $$\delta _M^{'}$$ = 0.35, nor with $$\delta _M^{'}$$ = 3.0. Fiddling with this parameter will certainly lead to an improved fit, and may allow the use of this approximative formula in certain cases, but the need to tune this parameter beforehand remains. 3.2 Case 2a: solid phase in rest, fluid phase moving In the next case, we assume that the fluid phase is moving with a constant velocity $$v_f^{'}$$ in z-direction. Due to non-dimensionalization, the velocity can directly be written as the Peclet number (i.e. $$v_f^{'} = \ {\rm{Pe}}$$) which can be defined by   \begin{equation}{\rm{Pe\ }} = \ \delta v \cdot {d_s} \cdot {\kappa ^{ - 1}},\end{equation} (30) with δv = vf − vs, i.e. in our case δv = vf. Since the z-direction is the only space coordinate, we can treat the problem in 1-D. We study the problem over a length scale of $$z_{{\rm{max}}}^{'} = \ 10$$. The initial conditions are $$T_f^{'}( {t{'} = \ 0} ) = T_s^{'}( {t{'} = \ 0} ) = \ 0$$ and the boundary conditions are $$T_f^{'}( {z{'} = \ 0} ) = T_s^{'}( {z{'} = \ 0} ) = \ 1$$ and $$T_f^{'}( {z{'} = z_{{\rm{max}}}^{'}} ) = T_s^{'}( {z{'} = z_{{\rm{max}}}^{'}} ) = \ 0$$. This condition implies a sudden heating at the bottom (z' = 0). The microscopic equivalent of this setup resembles Fig. 2(a) with hot fluid entering from below. Thus, our model describes the heating half-space problem including heat conduction and advection by porous flow. First, we explore the time evolution. Fig. 9(a) shows the solid and fluid temperatures assuming Pe = 3 and φ = 0.5 at times 0.25, 0.5, 1, 2 and 3. At the early times (0.25, 0.5), the pairs of curves resemble the erf function-like conductive cooling half-space solution, while at later times, advection takes over and the solutions develop into reverse s-shaped curves. This transitional behaviour stems from the fact that a conductively propagating temperature front moves with the velocity   \begin{equation}{v_{{\rm{diff}}}} = \ c\sqrt {\frac{\kappa }{t}} \end{equation} (31) with c of the order 1, while a convectively propagating front moves with the Darcy velocity $$\varphi v_f^{'}$$. Equating these velocities gives the overtake time, assuming c = 1)   \begin{equation}{t_{{\rm{ov}}}} = \frac{\kappa }{{v_f^2{\varphi ^2}}}\end{equation} (32) after which advection will take over and dominate the heat transport. This time can be written in non-dimensional form using the Peclet number to   \begin{equation}{t_{{\rm{ov}}}}^\prime = \frac{1}{{{\rm{P}}{{\rm{e}}^2}{\varphi ^2}}}\end{equation} (33) Figure 9. View largeDownload slide Temperature in the solid phase (blue lines) and the fluid phase (red lines) for successive times t'indicated in the figures, and for a setting with a porosity of 0.5. The fluid and the solid are suddenly heated up to T΄ = 1 at z΄ = 0 and the fluid moves upward in z’-direction with a velocity defined by the Peclet numbers (a) 3 and (b) 300. Figure 9. View largeDownload slide Temperature in the solid phase (blue lines) and the fluid phase (red lines) for successive times t'indicated in the figures, and for a setting with a porosity of 0.5. The fluid and the solid are suddenly heated up to T΄ = 1 at z΄ = 0 and the fluid moves upward in z’-direction with a velocity defined by the Peclet numbers (a) 3 and (b) 300. From this equation, it is clear that advection becomes important for the case shown in Fig. 9(a) for non-dimensional times larger than $${t_{{\rm{ov}}}}^{'} = \frac{4}{9}$$ (given Pe = 3 and φ = $$\frac{1}{2}$$), in agreement with the curve shapes changing from erf-type to S-type. Secondly, it is obvious, that from early times on the fluid is warmer than the solid. Thus, thermal non-equilibrium is clearly demonstrated. Hot fluid enters from below and is cooled by the solid, while the solid is heated by the advected fluid giving a temperature contrast along the curves of about 0.15 at early times decreasing slowly to 0.05 at t = 3. At later times, the temperature fronts of the fluid and solid move roughly with the same speed and the fluid temperature front is always about 0.25 diffusion lengths ahead. Thus, even on the macroscopic time and length scale such non-equilibrium temperatures are to be expected even for a Peclet number as low as 3. It is interesting to examine the evolution also in the limit of high Pe- numbers which is the relevant case for melt migrating through dykes. We carried out an experiment with Pe = 300 for φ = 0.5) and show the temperatures at the times 0.01, 0.02, 0.04, 0.08, 0.12, 0.16 and 0.2 in Fig. 9(b). One can clearly see the high-temperature front of the fluid entering and passing through the medium. Behind the front, the fluid has cooled only slightly to values above 0.9. The hot fluid acts as a continuous heating for the cold matrix within a region equivalent to the actual penetration depth of the fluid. Within this region, the matrix temperature continuously increases, while within the lowermost part (z' = 0–0.1) still a steep thermal boundary layer is present. Thus, a pronounced kink characterizes the resulting solid temperature profiles, separating the thin conductive thermal boundary layer and the deep region where penetrating fluid heats up the solid. Thus, for high Pe- numbers the maximum temperature difference between solid and fluid approaches 1. The lateral offset of the fluid and solid temperature curves increases indefinitely. In Fig. 10, we explore the effect of different Pe- numbers (fluid velocities) and porosities. The temperature in the fluid and the solid are shown at time t' = 1 (i.e. one diffusion time, see eq. 21) for an upward moving fluid for a porosity of φ = 0.1 (Fig. 10a) and φ = 0.5 (Fig. 10b), respectively. The results are shown for Pe- numbers equal to 0, 1, 3, 10 and 20 for φ = 0.5 (Fig. 10b) and for Pe- numbers equal 3, 10, 30 and 60 for φ = 0.1 (Fig. 10a). For the highest fluid velocity in each case (Pe = 60 for φ = 0.1 and Pe = 20 for φ = 0.5), the upper boundary condition at $$z_{{\rm{max}}}^{'}$$ strongly affects the results. We included the case for Pe = 0 (i.e. $$v_f^{'}$$ = 0), when the liquid is not moving, to verify that the solutions in the fluid and solid are identical and equal to the cooling (or heating) half-space solution (Turcotte & Schubert 2002) (as long as the heating front does not reach the right-hand side of the model). In fact, the analytical solution for the cooling half-space (which corresponds to Pe = 0) was also drawn in Fig. 10(b), and is perfectly matched by our solutions and thus not distinguishable. The time increments for the integration in eq. (26) depends on the porosity and were chosen according to the best solutions obtained for case 1 (see Fig. 8) to dt' = 10−5 for φ = 0.1 and to dt' = 10−4 for φ = 0.5. Figure 10. View largeDownload slide (a) and (b) Temperature in the fluid phase (red lines) and in the solid phase (blue lines) at the time t' = 1 (i.e. one diffusion time) for a setting where the fluid and the solid are suddenly heated up to T' = 1 at z' = 0 and the fluid is moving upward with different velocities (defined by the Peclet number Pe) for a porosity (a)φ = 0.1, (b)φ = 0.5, (c) non-dimensionalized temperature difference between fluid and solid phases and (d) spatial (i.e. vertical) offset of the fluid and solid temperature curves as a function of Peclet number times porosity (i.e. a measure for the volumetric flow rate). Note that the dimensionalized z is different in cases (a) and (b), because the scaling length ds is a function of φ. Figure 10. View largeDownload slide (a) and (b) Temperature in the fluid phase (red lines) and in the solid phase (blue lines) at the time t' = 1 (i.e. one diffusion time) for a setting where the fluid and the solid are suddenly heated up to T' = 1 at z' = 0 and the fluid is moving upward with different velocities (defined by the Peclet number Pe) for a porosity (a)φ = 0.1, (b)φ = 0.5, (c) non-dimensionalized temperature difference between fluid and solid phases and (d) spatial (i.e. vertical) offset of the fluid and solid temperature curves as a function of Peclet number times porosity (i.e. a measure for the volumetric flow rate). Note that the dimensionalized z is different in cases (a) and (b), because the scaling length ds is a function of φ. As can be expected, the temperature difference between the fluid and the solid phases increases with higher velocities (larger Peclet number) of the fluid. In case of a high volume fraction of the fluid (φ = 0.5, Fig. 10b), the heat content of the fluid is large and it heats the solid continuously up while moving upward. In case of a low volume fraction (φ = 0.1, Fig. 10a) its heat content is small and it is losing the heat quickly to the solid, rapidly cooling down, and the temperature in the large volume of the solid is only slightly increased. Thus higher Peclet numbers are needed to obtain a similar heating effect as for φ = 0.5. Interestingly, for similar penetration depths of the fluid and solid temperature fronts, their non-equilibrium temperature differences are similar and to first order independent of φ. The non-equilibrium temperature differences can be quantified by taking the temperature differences of the fluid and solid curves in Figs 10(a) and (b) at a given solid temperature (here 0.5) and plot the temperature difference (Fig. 10c) and the spatial distance between solid and fluid temperature values of 0.5 (Fig. 10d) as a function of the Darcy flow related Peclet number PeD  \begin{equation}{\rm{P}}{{\rm{e}}_D} = {\rm{Pe}} \cdot \varphi .\end{equation} (34) First, we observe that with increasing PeD the non-equilibrium temperature difference increases, but for the highest PeD the curves begin to flatten. It should be remembered that per definition the (non-dimensional) temperature difference is bounded by 1 and thus an indefinite increase is impossible (cf. Fig. 9b, which shows that for high Peclet number the temperature difference rapidly reaches values close to 1). Secondly, inspection of the time dependence of the temperature difference (e.g. shown in Fig. 9) indicates that the temperature difference decreases with time due to a decreasing temperature gradient with time. Thus, the curves shown in Fig. 10(c) show a snapshot, and for later times they will flatten. Thirdly, the spatial offset of the temperature curves (Fig. 10d) shows a clear linear relationship, depicted by a slope equal 1 in the log–log plot. Thus, for large PeD the spatial offset increases unboundedly. This is particularly important, as from inspection of the time behaviour (Fig. 9a), the offset of the curves does not show a detectable time-dependences, indicating that (i) even for times large compared to one diffusion time the spatial offset is still present, and (ii) the offset is approximately independent of the overall temperature gradient. 3.3 Comparison to an FE solution We also tested our formulation against a solution obtained from an explicit model of a propagating hot fluid through a set of parallel sheets in a colder surrounding. For this purpose, we created a finite-element model (FE, using COMSOL) representing a vertical section of non-dimensional length 10 of a fluid channel embedded in a resting solid. This setup represents case 2a on the local scale. Only half of one fluid layer in contact with half of a surrounding solid layer are modeled (with a resolution of 50 × 500), assuming symmetric thermal boundary conditions at both sides. The initial temperature was a set to zero everywhere. At the initial time, a temperature jump of amplitude 1 was prescribed at the lower boundary to both solid and fluid boundary sections. The fluid was assumed to move upwards with a constant non-dimensional velocity of Pe = 20 within the fluid layer. At the upper boundary, we used a zero-temperature gradient boundary condition both for the solid and fluid boundary sections. The thermal conductivities of fluid and solid are assumed to be equal. For this system, we solved the energy conservation equation, $$\frac{{\partial T{'}}}{{\partial t{'}}} + {\vec{v}_f}{'} \cdot \vec{\nabla }T{'}\ = \ {\nabla ^2}T{'}$$, with $${\vec{v}_f}{'}$$ abruptly changing to zero at the fluid–solid interface. The fluid fraction (i.e. porosity) was 50 per cent. The resulting temperature fields within the fluid and the solid phases are shown in Fig. 11(a) for different times. The hot fluid penetrates into the cold solid and generates a wavy- or finger-like hot temperature front which is diffusively broadening vertically with time, while the amplitude of the wave remains roughly constant. This is equivalent to the approximate time independence of the spatial distances of the fluid–solid temperature curves, for example, in Fig. 9(a). Figs 11(b)–(d) show the horizontally averaged temperatures within the solid and fluid phases in comparison with the temperatures determined by the macroscopic approach (eqs 24–26) using full integration. The agreement is very good. At the temperature front the numerical solutions seem to exhibit little numerical diffusion (see also Fig. 11a). At all times, the temperature differences between the numerical and macroscopic solutions are small, always less than 2 per cent. These differences result from both the resolution errors of the numerical FE model and the FD approximation of the heat eqs (23) and (25), as well as the numerical integration in eq. (26). Figure 11. View largeDownload slide (a) Three snapshots of the non-dimensional temperature field of an FE solution on the local scale. The computation domains have been mirrored four times for better visualization. (b)–(d) Comparison between the horizontally averaged temperatures (fluid and solid) obtained by the macroscopic approach as described in Section3.2 and by the FE model for a fluid sheet in a solid surrounding, resolving the fluid and the solid explicitly. Figure 11. View largeDownload slide (a) Three snapshots of the non-dimensional temperature field of an FE solution on the local scale. The computation domains have been mirrored four times for better visualization. (b)–(d) Comparison between the horizontally averaged temperatures (fluid and solid) obtained by the macroscopic approach as described in Section3.2 and by the FE model for a fluid sheet in a solid surrounding, resolving the fluid and the solid explicitly. 3.4 Comparison to approximate formula We also calculated the temperatures in both phases for a range of Peclet numbers (between 3 and 100) and fluid fractions φ (0.1 and 0.5) based on the approximate formula (eq. 27) and for different values of $$\delta _M^{'}$$ (i.e. the thickness of the microscopic thermal boundary layer) and compared the results to those obtained by the full backward in time integration. Fig. 12 shows the cases for (Pe, φ) equal to (30, 0.1) and (10, 0.5) and $$\delta _M^{'}$$ values of 0.05, 0.2, 0.3 and 0.5. Inspecting Fig. 12 suggests for the first case a $$\delta _M^{'}$$ value of about 0.2 and in the second case a value around 0.3. We calculated the rms deviations between the full integral solution and the approximate solutions for all cases and determined the minimum for the rms curves, which was $$\delta _M^{'}$$ (Pe = 30, φ = 0.1) = 0.2102 for the first case and $$\delta _M^{'}$$ (Pe = 10, φ = 0.5) = 0.3352 for the second case. As can be expected, in general we found a decrease in $$\delta _M^{'}$$ with higher Peclet numbers (in agreement with Amiri & Vafai 1994). However, we also found cases for small Peclet numbers where $$\delta _M^{'}$$ did not depend on Pe. Figure 12. View largeDownload slide Temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 based on the approximate eq. (27) (thick lines) and different values of $$\delta _M^{'}$$ in comparison with the results of the full backward integration in time (thin lines; same results as in Figs 10a and b) for fluid fraction φ = 0.1 and Peclet number Pe = 30 (upper row) and φ = 0.5 and Pe = 10 (lower row). Figure 12. View largeDownload slide Temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 based on the approximate eq. (27) (thick lines) and different values of $$\delta _M^{'}$$ in comparison with the results of the full backward integration in time (thin lines; same results as in Figs 10a and b) for fluid fraction φ = 0.1 and Peclet number Pe = 30 (upper row) and φ = 0.5 and Pe = 10 (lower row). To test whether the optimum $$\delta _M^{'}$$ changes with time (cf. Fig. 9 showing the time evolution of the full solution), we determined the optimum $$\delta _M^{'}$$ (Pe = 10, φ = 0.5) as 0.3352 at time t' = 1 and ranan approximate case with this $$\delta _M^{'}$$ for the time 0–3. The resulting evolution agrees extremely well (within line width) with the full solution (Fig. 9) for all times examined. This might be surprising because in case 1 we observed a significant disagreement between the approximate and the exact solution, however, this disagreement is most pronounced only for early times (t' < 0.25), while for later times as examined in Fig. 9, the instantaneous non-equilibrium approach seems to be sufficient if an appropriate $$\delta _M^{'}$$ is chosen. To summarize, the approximate formula (eq. 27) is a good approximation (and much simpler to implement in a numerical code for melt segregation) with a $$\delta _M^{'}$$value between 0.2 and 0.3 (depending on the average fluid (melt) velocity). 3.5 Case 2b: testing the necessary length of the convolution integral So far, we have carried out the integration in eq. (26) from the initial time $$t_0^{'}$$ = 0 to t' (i.e. along the entire fluid path). Inspecting the shape of the heat flow response function q* (Fig. 6) which rapidly decays for times of one diffusion time or less it might be expected that it is not necessary to carry out the convolution integral for the full time of evolution in order to save computation time. To check this, we defined a certain truncation time $$t_{{\rm{trunc}}}^{'} = t{'}\ - \Delta t_{{\rm{trunc}}}^{'}$$ and tested the quality of the solutions for the case φ = 0.1, a Peclet number of 3 and 10, and $$\Delta t_{{\rm{trunc}}}^{'}$$ = 0., 10−5, 10−4, 10−3, 10−2, 10−1 and 1. Results are plotted in Fig. 13 and show that the temperature in the fluid is strongly overestimated if the backward-in-time integration in eq. (26) is too short, while the temperature in the solid is much less affected. We conclude that always the full backward in time integrations is needed to avoid errors in solid–fluid heat transfer. The application to generalized two-phase flow requires a Lagrangian formulation to evaluate the integral and is thus computationally expensive. However, to avoid repeated calculations, it is advisable to store the result of the parentheses in the integrals in eqs (11) or (12) during the forward stepping. Figure 13. View largeDownload slide Dependency of the accuracy of the result on the length of the backward-in-time integration in eq. (26) for (a) slow fluid velocity (Pe = 3) and (b) faster fluid velocity (Pe = 10). Shown is the temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 and for a porosity of φ = 0.1. Figure 13. View largeDownload slide Dependency of the accuracy of the result on the length of the backward-in-time integration in eq. (26) for (a) slow fluid velocity (Pe = 3) and (b) faster fluid velocity (Pe = 10). Shown is the temperature in the fluid phase (red line) and in the solid phase (blue line) at the time t' = 1 (i.e. one diffusion time) for the same setting as in Fig. 10 and for a porosity of φ = 0.1. 4 DISCUSSION AND CONCLUSION 4.1 General remarks In this paper, we developed a theory to determine the heat exchange between a fluid and a solid phase on a macroscopic scale where the fluid is intruding or segregating into the solid and fluid and solid have initially different temperatures. The theory allows us to test the conditions under which different temperatures may prevail in the fluid and the solid phases. Such a formulation is a pre-requisite to model melt movement within solid rocks in regions where the temperature of the solid rock is below the solidus temperature. It may also be important when describing hydrothermal groundwater flow and heat transport through large fractures, which are not known well enough to be resolved as single structures. Heat exchange is based on a formula which implicitly assumes fluid-filled layers of equal spacing where the width of the layers depends on the volumetric fluid fraction. For a moving solid phase with a differentially moving fluid phase we derived an expression to calculate the heat exchange rate across the fluid–solid interfaces using this formula and integrating the heat exchange history along solid flow path. Solving this integral gives an interphase heat exchange rate of dimension Jm−3 s−1, which couples the two heat transport equations for the solid and the fluid phases. Solving these heat equations allows determining the spatially and temporarily varying temperatures in the fluid and the solid phases separately. 4.2 Testing thermal non-equilibrium by using Peclet numbers Before discussing the possible applications of our formulation, it is worthwhile to introduce alternative Peclet numbers which may be more appropriate for various scenarios. While our Peclet number (eq. 30) is based on the fluid velocity and a length scale describing the distance ds of fluid inclusions (dykes and pores), a macroscopic Peclet number based on the length scale of the system L may be defined by   \begin{equation}{\rm{P}}{{\rm{e}}_L} = \frac{{\delta vL}}{\kappa }\quad {\rm{with}}\quad {\rm P}{{\rm e}_L} = \ {\rm{Pe}}\frac{L}{{{d_s}}}\end{equation} (35) Often the volumetric flow rate (Darcy velocity vD) is used instead of the fluid velocity leading to the Darcy related Peclet number   \begin{equation}{\rm{P}}{{\rm{e}}_D} = \frac{{{v_D}{d_s}}}{\kappa }\quad {\rm{with}}\quad {\rm{P}}{{\rm{e}}_D} = \ {\rm{Pe}}\ \varphi \quad {\rm{from}}\,{\rm{eq}}.\,\,(34).\end{equation} (36) It is noteworthy that the overtaking time tov' at which Darcy related advection takes over with respect to heat conduction (cf. eq. 33) is directly given by PeD− 2. If the Peclet number is based on the Darcy velocity and a macroscopic length scale of the system L, we obtain   \begin{equation}{\rm{P}}{{\rm{e}}_{LD}} = \frac{{{v_D}L}}{\kappa }\quad {\rm{with}}\quad {\rm{P}}{{\rm{e}}_{LD}} = \ {\rm{Pe}}\ \frac{{\varphi L}}{{{d_s}}}\end{equation} (37) In many applications of two-phase flow, neither the microscopic, nor the macroscopic length scales are used, but the problems are scaled by the compaction length (McKenzie 1984). Using McKenzie's definition, the compaction length δc can be written in terms of our characteristic length scale ds by   \begin{equation} {\delta _c} = {\left( {\frac{{\left( {{\eta _b} + \frac{4}{3}{\eta _s}} \right){\varphi ^n}}}{{{\eta _f}b}}} \right)^{{\raise0.7ex\hbox{1} / \lower0.7ex \hbox{2}}}}\ {d_s}. \end{equation} (38) Here ηb, ηs and ηf are the bulk, shear and fluid viscosities, respectively, b is a geometric constant of the order 100 and n is the power of the permeability–porosity relation. For realistic partially molten systems, the compaction length is thus 2–7 orders of magnitude larger than ds. The main results of this paper may be used to check for a given Darcy flow scenario whether the temperature field is dominated by Darcy advection and whether fluid and solid temperatures are different due to thermal non-equilibrium. A quick check starts from determining the characteristic Peclet number and porosity of the system to obtain PeD. If the characteristic non-dimensional timescale of the flow problem is larger than PeD− 2, then the temperature field will be dominated by Darcy flow rather than conduction alone. Furthermore, Figs 10(c) and (d) imply that in the presence of an overall temperature gradient in fluid flow direction fluid and solid are not in thermal equilibrium if PeD > O(1) to O(10). Under such conditions, the temperatures may be different by 10’s per cent of the typical temperature variations of the problem. This is equivalent to the fact that the same temperatures in the solid and fluid phases may be spatially separated by a distance which is proportional to the characteristic length of the solid phase, ds, times PeD. 4.3 Combining the thermal non-equilibrium formulation with general two-phase flow In the test cases as presented in this paper, the fluid velocity is prescribed. However, in a hydrothermal flow problem with a given rock permeability distribution, the fluid flow can be calculated by solving additionally the Darcy and mass conservation equations. The convolution integral over time (eq. 26) still needs to be evaluated for every point in space. In the full two-phase flow case, the matrix may be deformable and may flow according to solutions of the Stokes equations. To apply our thermal non-equilibrium formulation first, the flow paths of the matrix (not of the fluid) need to be determined and stored. At a given time step, the convolution integral must be carried out backward in time along the matrix path, which implies that the temperature changes need to be known everywhere for all times before the given time step. While rigorous thermal non-equilibrium modeling of porous flow requires the numerically expensive full backwards integrations of the convolution integral (which means that the temperature fields for solid and fluid phases have to be stored for each elapsed time step and considered in the convolution integral for each new time step), we showed that for times larger than the characteristic diffusion time the approximate formula (eq. 27) (cf. approximations by Amiri & Vafai 1994 or Minkowycz et al. 1999) does a good job and significantly reduces computation times. This is particularly important for full two-phase flow formulations in which our theory would imply that also the evolving flow paths of the matrix need to be tracked and integrated backward in time. Thus, if this integration could be replaced by an approximate instantaneous heat exchange term, two-phase flow formulations such as shown by, for example, Schmeling (2000) or Wallner and Schmeling (2016) could be easily extended by separating the heat equation into two equations and coupling them by the instantaneous heat exchange term. Our method may also be applied to fluid geometries others than regularly spaced channels or dykes. This requires alternative formulations of the heat flow response function to a temperature step function, $$q_s^*( t ),$$ that is, replacement of eq. (17). For more complicated geometries, such functions can be determined numerically by solving an initial value problem of heat conduction and storing $$q_s^*( t )$$ in a lookup table. We also want to emphasize that thermal non-equilibrium in hydrothermal porous flow problems is generally of transient nature, arising from a sudden change in the physical state: for example, a sudden change in the flow as when a dyke is injected. Thus, our formulation needs to be lined to a formalism which calculates of estimates the changing flow field. In case of steady state or slowly changing flow fields, sudden changes in the physical state may include porous-driven temperature fronts, in which thermal non-equilibrium is predicted to occur in the high-temperature gradient regions as shown in this paper. 4.4 Application to melt segregation at the transition between channeling and dykes One of the most obvious two-phase flow systems where our theory may be applied is the transitional regime between melt channeling within a melting source zone (Stevenson 1989; Aharonov et al. 1995) and magma ascent by dykes through the subsolidus mantle lithosphere or crust (e.g. Rivalta et al.2015). Keller et al. (2013) studied isothermal melt ascent by numerical models, focusing on rheology but neglected freezing, and identified regimes of viscous diapirism, viscoplastic decompaction channels and elastoplastic dyking. Leitch and Weinberg (2002) modeled mesoscale pervasive magma ascent through dykes but assumed thermal equilibrium between dykes and country rocks. Our results allow testing this kind of assumption. We will make a short, incomplete tour through these processes with the emphasis of extracting the parameters controlling the Darcy related Peclet number PeD, that is, we need estimates about the channel or dyke distance, ds and the volumetric flow rate vD = vf φ. See also Fig. 14 for illustration of this tour. Figure 14. View largeDownload slide Illustration of the discussion of melt ascent and associated Peclet numbers which result in predictions about thermal equilibrium or non-equilibrium during this ascent process. The horizontal solidus temperature indicates the situation, where melt temperatures will be slightly below solidus temperatures, thus the melt has the potential to penetrate into the subsolidus region. Figure 14. View largeDownload slide Illustration of the discussion of melt ascent and associated Peclet numbers which result in predictions about thermal equilibrium or non-equilibrium during this ascent process. The horizontal solidus temperature indicates the situation, where melt temperatures will be slightly below solidus temperatures, thus the melt has the potential to penetrate into the subsolidus region. From the analysis of Stevenson (1989), a partially molten rock undergoing deviatoric straining develops a channeling instability with a characteristic wavelength equal or smaller (down to several grain sizes) than the compaction length, which is defined by $${( {{k_\varphi }( {{\eta _b} + \frac{{4{\eta _s}}}{3}} )/{\eta _f}} )^{1/2}}$$, with kφ as initial permeability, ηb,  and  ηs as bulk and shear viscosity of the matrix and ηf as melt viscosity. Typical values for partially molten mantle give compaction lengths between 10s of metres and few kilometres. Models on channel formation by Richardson (1998) and Golabek et al. (2008) found channeling instabilities on larger scales, but these were probably limited to the finite grid resolution, inhibiting the development of smaller wavelengths. Well-resolved models of vein formation due to shearing and compaction by Rabinowicz & Toplis (2009) found the formation of regular veins with distances of 0.03 compaction lengths (several 10s of metres). Modeling of channeling instabilities due to reaction infiltration in compacting media by Spiegelman et al. (2001) found initial distances between veins of about 0.3 compaction lengths, but upon ascending they continuously merge into fewer channels with larger distances. Such coalescence of melt conduits has been proposed as the ‘fractal tree model’ by Hart (1993) based on geochemical arguments. Observation (e.g. Weinberg et al. 2015) shows that the distance between veins may be on the order of decimetres. In summary, we expect that the channel distance close to the onset of channeling within the melt source regions may be rather small, on the order of decimetres, but further up, still in the supersolidus partially molten regions, may increase up to the metre or kilometre scale. The segregation velocity of magma in highly permeable or fully melt-filled veins depends mainly on the melt viscosity. A velocity of around 2.5 m yr−1 (Petford 1995) has been estimated as an average value, however, this will strongly depend on the spacing and widths of the channels. We might therefore better try to estimate the volumetric flow rates (Darcy velocity). Let us take as a typical example of melt ascent a mid-ocean ridge, where a basaltic crust of thickness h spreads with a half-spreading rate of v0. If the melt is collected within an asthenospheric wedge of half-width L (which, as a consequence of the wedge shape, decreases towards the surface), on first order the Darcy velocity within the mantle wedge is given by v0h/L. For typical mid-ocean ridges, this results in a volumetric melt flow rate of 1e-10 to 1e-9 ms−1, the smaller value being more typical for greater depth (where L is wider). Such values have also been reported from the fully dynamic melting models by Keller et al. (2017). With a thermal diffusivity κ = 10−6 m2 s−1, a small (because deep) Darcy velocity and dense veins, the resulting Pe- numbers (PeD) are of the order of 10−5–10−3, that is, the flow is in thermal equilibrium. At shallower depth, and when melt conduits have merged, PeD may increase to order 1 or larger, that is, in the presence of a vertical temperature gradient thermal non-equilibrium effects will become moderately important (see Fig. 14). It should be noted that these rough estimates are conservative estimates, because they assume (1) steady-state ascent of melt and (2) passive upwelling. Regions with active mantle upwelling and time-dependent melt extraction will exhibit higher Peclet numbers implying a high potential for thermal non-equilibrium. When the above-mentioned melt channels reach the level of the solidus temperature, they must undergo a transition into dykes. Dykes are generally believed to be generated by magmatic overpressure and brittle failure of rocks and propagate through initially unfractured rock (Rivalta et al. 2015). We will now focus on this end-member of magma ascent and estimate appropriate Peclet numbers. The sizes of dykes in nature is very different, however, a typical dyke has a thickness between some 10s of centimetres to some 10s of metres with a length of a few 100 m. In the vertical plane, dykes (or in the horizontal plane for sills) can be some tens of kilometres (Wright et al. 2006). There is also a large amount of literature reporting the occurrence of tens to thousands of subparallel dykes and/or melt veins. Such dykes can be found, for example, in the vicinity of active volcanoes (Gudmundsson 1983; 2002), in ophiolites formed near mid-oceanic ridges (e.g. Nicolas et al. 1994) or in continental volcanic provinces (e.g. Ray et al. 2007) and the magmatic sources are proposed in the upper or lower crust (e.g. Guernina & Sawyer 2003) or even in the upper mantle (Ernst et al. 2001). Whether multiple dykes in a swarm are propagating at the same time is still a matter of investigation (e.g. Jin & Johnson 2008). Havlin et al. (2013) modeled dyke initiation from a melt accumulation layer at the top of the asthenosphere and obtained dyke recurrence times of the order of days and spacings of 100 m or larger. Field observations show at least that magma propagation in adjacent dykes during active rifting events occur within short time intervals of a few months (e.g. Einarsson & Brandsdottir 1980; Grandin et al. 2010). Dahm (2000) estimated for an average dyke with a width of a few decimetre and a length of a few kilometres a propagation velocity (which is equivalent with the propagation of the hot magmatic fluid) on the order of 0.1–1 ms−1: recalling the definition of the Peclet number (eq. 30) and assuming a dyke spacing of 100 m leads to a Peclet number on the order of 105, thus by far larger than studied in this papers, indicating a strong non-equilibrium temperature distribution between fluid and solid. From this exercise, we conclude (Fig. 14) that within the source regions of mantle melting steady-state porous flow through a channelized system is in thermal equilibrium with the matrix. Once this flow becomes episodic, with volumetric flow rates exceeding the steady-state rate be one or two orders of magnitude, the Peclet number becomes so large that thermal non-equilibrium effects may become important. Furthermore, within the transition zone between melt rising through highly permeable channels and melt rising through propagating dykes the Peclet number is expected to strongly exceed the value one. Within this transition zone, our results indicate that melt and matrix temperature may be different by an order of 10 per cent of the characteristic temperature variations of the tectonic setup, implying that melt is still above the solidus, while the solid is below. Entering the dyke regime, extremely large Peclet numbers are predicted implying large temperature differences between rock and fluid. 4.5 Summary We conclude with summarizing the main key points of this study: A formulation of calculating thermally non-equilibrium two-phase flow has been developed based on the time-dependent heat exchange between fluid and solid which integrates the thermal history along flow paths by a convolution integral. For Darcy velocity-based Peclet numbers larger 1 thermal non-equilibrium between fluid and matrix becomes important with temperature differences of 5 per cent or more relative to characteristic temperature differences of the flow problem. The numerically expensive calculation of the history-dependent heat exchange can be replaced in the long-term limit by using an approximate instantaneous exchange term with adequate parameters. Our results imply that thermal non-equilibrium can play an important role for melt migration through partially molten systems where melt focuses into melt channels and near the transition to melt ascent by dykes. Acknowledgements This paper was accomplished during a sabbatical stay of the authors at the Monash University, Melbourne. We are grateful to all our colleagues from the School of Earth, Atmosphere and Environment for their hospitality and support. We thank two anonymous reviewers for their patience to go through all equations and to contribute several constructive comments, which helped us to improve the paper considerably. We acknowledge funding support for this sabbatical stay by the Deutsche Forschungsgemeinschaft (DFG) with the grant no. SCHM 872/25-1. REFERENCES Aharonov E., Whitehead J.A., Kelemen P.B., Spiegelman M., 1995. Channeling instability of upwelling melt in the mantle. J. geophys. Res. , 100( 20), 433−450. Amiri A., Vafai K., 1994. Analysis of dispersion effects and non-thermal equilibrium, non-Darcian, variable porosity incompressible flow through porous media, Int. J. Heat Mass Transf. , 37 939– 954. Google Scholar CrossRef Search ADS   Bercovici D., Ricard Y., Schubert G., 2001. A two-phase model for compaction and damage, Part 1: general theory, J. geophys. Res. , 106( B5), 8887– 8906. Google Scholar CrossRef Search ADS   Dahm T., 2000. On the shape and velocity of fluid-filled fractures in the Earth, Geophys. J. Int . 142 181– 192. Google Scholar CrossRef Search ADS   de Lemos M.J.S., 2016. Thermal Non-equilibrium in Heterogeneous Media , Springer Science+Business Media, Inc. Google Scholar CrossRef Search ADS   Golabek G.J., Schmeling H., Tackley P.J., 2008. Earth's core formation aided by flow channeling instabilities induced by iron diapirs, Earth planet. Sci. Lett.   271 24– 33. Google Scholar CrossRef Search ADS   Grandin R., Socquet A., Jacques E., Mazzoni N., de Chebalier J.B., King G., 2010. Sequence of rifting in Afar, Manda-Hararo rift, Ethiopia, 2005–2009: time–space evolution and interactions between dikes from interferometric synthetic aperture radar and static stress change modeling, J. geophys. Res. , 115 B10413, doi:10.1029/2009JB000815. Google Scholar CrossRef Search ADS   Gudmundsson A., 1983. Form and dimensions of dykes in eastern Iceland, Tectonophysics , 95 295– 307. Google Scholar CrossRef Search ADS   Gudmundsson A., 2002. Emplacement and arrest of sheets and dykes in central volcanoes, J. Volc. Geotherm. Res. , 116 279– 298. Google Scholar CrossRef Search ADS   Guernina S., Sawyer E.W., 2003. Large-scale melt-depletion in granulite terranes: an example from the Archean Ashuanipi Subprovinces of Quebec, J. Metamorphic Geol. , 21 181– 201. Google Scholar CrossRef Search ADS   Einarsson P., Brandsdottir B. 1980. Seismological evidence for lateral magma intrusion during the July 1978 deflation of the Krafla volcano in NE Iceland, J. Geophys. , 47 160– 165. Ernst R.E., Grosfils E.B., Mège D., 2001. Giant dike swarms: Earth, Venus, and Mars. Annu. Rev. Earth planet. Sci. , 29 489– 534. Google Scholar CrossRef Search ADS   Hart S.R., 1993. Equilibrium during melting: a fractal tree model, Proc. Natl. Acad. Sci. USA , 90 11914– 11918. Google Scholar CrossRef Search ADS   Havlin C., Parmentier E.M., Hirth G., 2013. Dike propagation driven by melt accumulation at the lithosphere asthenosphere boundary, Earth planet. Sci. Lett. , 376 20– 28. Google Scholar CrossRef Search ADS   Jin Z.-H., Johnson S.E.. 2008. Magma-driven multiple dike propagation and fracture toughness of crustal rocks, J. geophys. Res. , 113 B03206, doi:10.1029/2006JB004761. Google Scholar CrossRef Search ADS   Katz R.F., 2008. Magma dynamics with the enthalpy method: benchmark solutions and magmatic focusing at mid‐ocean ridges, J. Petrol. , 49( 12), 2099– 2121. Google Scholar CrossRef Search ADS   Keller T., May D.A., Kaus B.J.P., 2013. Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust, Geophys. J. Int. , 195( 3), 1406– 1442. Google Scholar CrossRef Search ADS   Keller T., Katz R.F., Hirschmann M.M., 2017. Volatiles beneath mid-ocean ridges: deep melting, channelised transport, focusing, and metasomatism, Earth planet. Sci. Lett. , 464 55– 68. Google Scholar CrossRef Search ADS   Leitch A.M., Weinberg R.F., 2002. Modellinggranite migration by mesoscale pervasive flow, Earth planet. Sci. Lett. , 200 131– 146. Google Scholar CrossRef Search ADS   MATLAB and Statistics Toolbox Release 2012b, The MathWorks, Inc., Natick, MA. McKenzie D., 1984. The generation and compaction of partially molten rock, J. Petrol. , 25 713– 765. Google Scholar CrossRef Search ADS   Minkowycz W.J., Haji-Sheikh A., Vafai K., 1999. On departure from local thermal equilibrium in porous media due to a rapidly changing heat source: the Sparrow number, Int. J. Heat Mass Transf. , 42 ( 18), 3373– 3385. Google Scholar CrossRef Search ADS   Nield D.A., Bejan A., 2006. Convection in Porous Media , 3rd edn, Springer Science+Business Media, Inc. Nicolas A., Boudier F., Ildefonse B., 1994. Evidence from the Oman ophiolite for active mantle upwelling beneath a fast spreading ridge, Nature , 370 51– 53. Google Scholar CrossRef Search ADS   Petford N., 1995. Segregation of tonalitic-trondhjemitic melts in the continental crust: the mantle connection, J. geophys. Res. , 100 15735– 15743. Google Scholar CrossRef Search ADS   Rabinowicz M., Toplis M.J., 2009. Melt segregation in the lower part of the partially molten mantle zone beneath an Oceanic Spreading Centre: numerical modelling of the combined effects of shear segregation and compaction, J. Petrol. , 50 ( 6), 1071– 1106. Google Scholar CrossRef Search ADS   Ray R., Sheth H.C., Mallik J., 2007. Structure and emplacement of the Nandurbar–Dhule mafic dyke swarm, Deccan Traps, and the tectonomagmatic evolution of flood basalts, Bull. Volc ., 69 537– 551. Google Scholar CrossRef Search ADS   Richardson C.N., 1998. Melt flow in a variable viscosity matrix, Geophys. Res. Lett. , 25 1099– 1102. Google Scholar CrossRef Search ADS   Rivalta E., Taisne B., Buger A.P., Katz R.F. 2015. A review of mechanical models of dike propagation: schools of thought, results and future directions, Tectonophysics , 638 1– 42. Google Scholar CrossRef Search ADS   Rubin A.M., 1998. Dike ascent in partially molten rock, J. geophys. Res. , 103 ( B9), 20901– 20919. Google Scholar CrossRef Search ADS   Rudge J.F., Bercovici D., Spiegelman M., 2011. Disequilibrium melting of a two phase multicomponent mantle, Geophys. J. Int. , 184( 2), 699– 718. Google Scholar CrossRef Search ADS   Schmeling H., 2000. Partial melting and melt segregation in a convecting mantle, in Physics and Chemistry of Partially Molten Rocks , eds. Bagdassarov N., Laporte D., Thompson A.B., pp. 141– 178. Kluwer Academic Publ., Dordrecht. Google Scholar CrossRef Search ADS   Spiegelman M., 1993. Flow in deformable porous media. Part 1. Simple analysis, J. Fluid Mech. , 247 17– 38. Google Scholar CrossRef Search ADS   Spiegelmann M., Kelemen P.B., Aharonov, E., 2001. Causes and consequences of flow organization during melt transport: the reaction infiltration instability in compactible media, J. geophys. Res. , 106 2061– 2077. Google Scholar CrossRef Search ADS   Stevenson D.J., 1989. Spontaneous small-scale melt segregation in partial melts undergoing deformation, Geophys. Res. Lett. , 16( 9), 1067– 1070. Google Scholar CrossRef Search ADS   Turcotte D.L., Schubert G., 2002. Geodynamics . Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Wallner H., Schmeling H., 2016. Sensitivity analysis of rift induced delamination with application to Rwenzori Mountains, Geophys. J. Int. , 187 1135– 1145. Google Scholar CrossRef Search ADS   Weinberg R.F., Veveakis E., Regenauer-Lieb K., 2015. Compaction-driven melt segregation in migmatites, Geology , doi:10.1130/G36562.1. Wright T., Ebinger C., Biggs J., Ayele A., Yirgu G., 2006. Magma-maintained rift segmentation at continental rupture in the 2005 Afar dyking episode, Nature , 442 291– 294. Google Scholar CrossRef Search ADS PubMed  © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off