2.3.1 Local energy conservation

In the laminar regime and in axisymmetrical coordinates ( $r$ , $z$ ), the local energy conservation equation within the lubrication assumption is

(2.7) $$\begin{eqnarray}\displaystyle \frac{\text{D}}{\text{D}t}(\unicode[STIX]{x1D70C}_{m}C_{p,m}T+\unicode[STIX]{x1D70C}_{m}L(1-\unicode[STIX]{x1D719}))=k_{m}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}, & & \displaystyle\end{eqnarray}$$

where $T(r,z,t)$ is the fluid temperature and $\unicode[STIX]{x1D70C}_{m}$ , $k_{m}$ and $C_{p,m}$ are the density, thermal conductivity and specific heat of the fluid. Here, we also account for energy release by possible crystallization of the fluid, which is a non-negligible source of heat in the case of magmas; $\unicode[STIX]{x1D719}(r,z,t)$ is the crystal fraction in the melt and $L$ the latent heat of crystallization. In this model, the crystals are only considered as a source/sink of energy as they melt/form at equilibrium during the flow as we assume that the physical properties of the crystal liquid mixture are the same as that of the fluid.

The fluid temperature varies between its liquidus $T_{L}$ and solidus temperature $T_{S}$ , i.e. $T_{i}=T_{L}$ and $T_{0}=T_{S}$ . As $T_{0}$ is also the fixed temperature of the surrounding medium in this model, this is equivalent to considering a fluid/wall interface pinned near the bulk freezing temperature $T_{S}$ of the fluid. In making this assumption, we neglect the heat flux necessary to heat up the wall up to $T_{S}$ .

Following a common approximation, we assume that the crystal fraction is a linear function of temperature over the melting interval (Hort Reference Hort1997; Michaut & Jaupart Reference Michaut and Jaupart2006)

(2.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=\frac{T_{L}-T}{T_{L}-T_{s}}=\frac{T_{i}-T}{T_{i}-T_{0}}. & & \displaystyle\end{eqnarray}$$

With these approximations, the local energy equation (2.7) becomes

(2.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}r}+w\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}=\frac{St}{St+1}\unicode[STIX]{x1D705}_{m}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}, & & \displaystyle\end{eqnarray}$$

where $u(r,z,t)$ and $w(r,z,t)$ are the radial and vertical fluid velocities, $St=(C_{p,m}(T_{i}-T_{0}))/L$ is the Stefan number which represents the ratio of sensible heat between solidus and liquidus to the total energy of the fluid at the liquidus temperature and $\unicode[STIX]{x1D705}_{m}$ is the fluid thermal diffusivity $\unicode[STIX]{x1D705}_{m}=k_{m}/(\unicode[STIX]{x1D70C}_{m}C_{p,m})$ . In particular, equation (2.9) shows that considering the energy released by crystallization is simply equivalent to considering a reduced thermal diffusivity $\widetilde{\unicode[STIX]{x1D705}_{m}}$ , which reads

(2.10) $$\begin{eqnarray}\displaystyle \widetilde{\unicode[STIX]{x1D705}_{m}}=\frac{St}{St+1}\unicode[STIX]{x1D705}_{m}, & & \displaystyle\end{eqnarray}$$

in the heat transport equation.

We use an integral balance method of heat transfer theory to approximately solve equation (2.9). In this method, that is developed below, the vertical structure of the temperature field is represented by a known function of depth that approximates the expected solution (Goodman Reference Goodman1958). Previous works on the cooling of lava domes at the surface have shown that such a reduction efficiently reduces the computation time while keeping the full dynamics of the unsimplified equation (2.9) well resolved (Balmforth et al. Reference Balmforth, Craster and Sassi2004). We use different kinds of vertical temperature profiles and show in appendix A that they all lead to very similar results.

2.3.2 Integral balance solution for the temperature $T(r,z,t)$

We model the cooling of the flow through the growth of two thermal boundary layers: one growing downward from the top and a second growing upward from the base. As we assume a fixed temperature at the boundary, the two thermal boundary layers grow symmetrically and have the same thickness $\unicode[STIX]{x1D6FF}(r,t)$ (figure 1). A popular approximation for the vertical temperature profile $T(r,z,t)$ is

(2.11) $$\begin{eqnarray}\displaystyle T=\left\{\begin{array}{@{}ll@{}}T_{b}-(T_{b}-T_{0})\left(1-\displaystyle \frac{z}{\unicode[STIX]{x1D6FF}}\right)^{n}\quad & 0\leqslant z\leqslant \unicode[STIX]{x1D6FF}\\ T_{b}\quad & \unicode[STIX]{x1D6FF}\leqslant z\leqslant h-\unicode[STIX]{x1D6FF}\\ T_{b}-(T_{b}-T_{0})\left(1-\displaystyle \frac{h-z}{\unicode[STIX]{x1D6FF}}\right)^{n}\quad & h-\unicode[STIX]{x1D6FF}\leqslant z\leqslant h,\end{array}\right. & & \displaystyle\end{eqnarray}$$

where $T_{b}(r,t)$ is the temperature at the centre of the flow and $n>1$ (Balmforth et al. Reference Balmforth, Craster and Sassi2004). This approximation captures the essential behaviour of the thermal structure: cooling is concentrated at the upper and bottom interfaces and cold boundary layers grow into the fluid interior as it flows. In addition, this profile assures the continuity of the temperature and heat flux within the flow.

While higher-order contributions to the temperature field, i.e. with $n>2$ , may also exist in certain situations, a parabolic profile is the most natural choice and we use $n=2$ in the following (Bercovici Reference Bercovici1994; Bercovici & Lin Reference Bercovici and Lin1996). Nonetheless, the case of higher-order profiles are treated in appendix A and are shown not to influence the flow dynamics at least within the level of description adopted here.

While the integral balance solution in (2.11) depends on two variables $T_{b}$ and $\unicode[STIX]{x1D6FF}$ that have to be consistently determined, these two variables are not independent. Indeed, in the flow region where thermal boundary layers exist, i.e. where $\unicode[STIX]{x1D6FF}<h/2$ , the temperature $T_{b}$ is equal to the injection temperature $T_{i}$ and $\unicode[STIX]{x1D6FF}$ is the unknown (figure 1, $r=r_{b}$ ). In contrast, in the flow region where the thermal boundary layers have connected, which eventually arises at the current front, $\unicode[STIX]{x1D6FF}=h/2$ and $T_{b}$ becomes the unknown (figure 1, $r=r_{a}$ ).

2.3.3 Integral balance equation

We integrate the local energy conservation equation (2.9) separately over the two thermal boundary layers. The integration over the bottom thermal layer, from $z=0$ to $z=\unicode[STIX]{x1D6FF}$ gives

(2.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D6FF}(\bar{T}-T_{b}))+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\unicode[STIX]{x1D6FF}(\overline{uT}-\bar{u}T_{b}))+\unicode[STIX]{x1D6FF}\left(\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}t}+\overline{u}\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =-\widetilde{\unicode[STIX]{x1D705}_{m}}\left.\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}\right|_{z=0}+w_{i}(T_{i}-T_{b}),\end{eqnarray}$$

where the bar indicates a vertical average over the bottom thermal boundary layer

(2.13) $$\begin{eqnarray}\displaystyle \overline{f}=\frac{1}{\unicode[STIX]{x1D6FF}}\int _{0}^{\unicode[STIX]{x1D6FF}}f\,\text{d}z, & & \displaystyle\end{eqnarray}$$

$T_{b}(r,t)$ is the temperature at $z=\unicode[STIX]{x1D6FF}$ and we have used the nullity of the thermal gradient at $z=\unicode[STIX]{x1D6FF}$ and the local mass conservation

(2.14) $$\begin{eqnarray}\displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}ru}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}=0. & & \displaystyle\end{eqnarray}$$

The integration over the top thermal layer, from $z=h-\unicode[STIX]{x1D6FF}$ to $z=h$ gives

(2.15) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D6FF}(\bar{T}-T_{b}))+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\unicode[STIX]{x1D6FF}(\overline{uT}-\bar{u}T_{b}))+\unicode[STIX]{x1D6FF}\left(\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}t}+\overline{u}\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\widetilde{\unicode[STIX]{x1D705}_{m}}\left.\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}\right|_{z=h},\end{eqnarray}$$

where, in addition to (2.14) and the fact that the thermal gradient at $z=h-\unicode[STIX]{x1D6FF}$ is equal to zero, we have used the kinematic boundary condition at $z=h(r,t)$

(2.16) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}=w. & & \displaystyle\end{eqnarray}$$

Adding (2.12) and (2.15) and using (2.11) with $n=2$ to derive the conductive fluxes, we finally obtain the heat balance equation which reads

(2.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D6FF}(\bar{T}-T_{b}))+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\unicode[STIX]{x1D6FF}(\overline{uT}-\bar{u}T_{b}))+\unicode[STIX]{x1D6FF}\left(\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}t}+\overline{u}\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =-2\widetilde{\unicode[STIX]{x1D705}_{m}}\frac{(T_{b}-T_{0})}{\unicode[STIX]{x1D6FF}}+\frac{w_{i}}{2}(T_{i}-T_{b}).\end{eqnarray}$$

2.3.4 Final heat transport equation

Equation (2.17) can be rewritten depending on the unknown in the temperature field $\unicode[STIX]{x1D6FF}$ or $T_{b}$ .

When $T_{b}=T_{i}$ and $\unicode[STIX]{x1D6FF}$ is the variable, equation (2.17) becomes

(2.18) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D6FF}(\bar{T}-T_{i}))+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\unicode[STIX]{x1D6FF}(\overline{uT}-\bar{u}T_{i}))=-2\widetilde{\unicode[STIX]{x1D705}_{m}}\frac{(T_{i}-T_{0})}{\unicode[STIX]{x1D6FF}}. & & \displaystyle\end{eqnarray}$$

When $\unicode[STIX]{x1D6FF}=h/2$ and $T_{b}$ is the variable, equation (2.17) becomes

(2.19) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h(\bar{T}-T_{b}))+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(rh(\overline{uT}-\bar{u}T_{b}))+h\left(\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}t}+\overline{u}\frac{\unicode[STIX]{x2202}T_{b}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =-8\widetilde{\unicode[STIX]{x1D705}_{m}}\frac{(T_{b}-T_{0})}{h}+w_{i}(T_{i}-T_{b}).\end{eqnarray}$$

Expanding the derivatives and using a global statement of mass conservation, equation (2.19) simplifies and reads

(2.20) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h\bar{T})+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(rh\overline{uT})=-8\widetilde{\unicode[STIX]{x1D705}_{m}}\frac{(T_{b}-T_{0})}{h}+w_{i}T_{i}. & & \displaystyle\end{eqnarray}$$

Interestingly, equation (2.20) is a particular case of (2.18) if the boundary layers have merged and $\unicode[STIX]{x1D6FF}=h/2$ . We can then use (2.18) as a simplification of (2.17) to describe the heat transport within the flow in both cases.

We finally rewrite (2.18) using a new variable $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FF}(T_{i}-\bar{T})$ which encloses both unknowns $\unicode[STIX]{x1D6FF}$ and $T_{b}$ since $\bar{T}$ depends on $T_{b}$ following $\bar{T}=(2T_{b}+T_{0})/3$

(2.21) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}t}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\bar{u}\unicode[STIX]{x1D709})-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\unicode[STIX]{x1D6FF}(\overline{uT}-\bar{u}\bar{T}))=2\widetilde{\unicode[STIX]{x1D705}_{m}}\frac{(T_{i}-T_{0})}{\unicode[STIX]{x1D6FF}}. & & \displaystyle\end{eqnarray}$$

The second term on the left-hand side of (2.21) contains advection by the vertically integrated radial velocity profile while the third term contains a correction accounting for the vertical structure of the temperature field. The term on the right is the loss of heat by conduction in the surrounding medium.

In addition, this formulation in terms of a unique variable $\unicode[STIX]{x1D709}$ allows us to calculate $T_{b}$ or $\unicode[STIX]{x1D6FF}$ directly from the expression of $\unicode[STIX]{x1D709}$ using either $\unicode[STIX]{x1D6FF}=h/2$ or $T_{b}=T_{i}$ respectively

(2.22a,b ) $$\begin{eqnarray}\displaystyle T_{b}(r)=\left\{\begin{array}{@{}ll@{}}T_{i}\quad & \text{if }\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x1D709}_{t}\\ \displaystyle \frac{3T_{i}-T_{0}}{2}-\frac{3\unicode[STIX]{x1D709}}{h}\quad & \text{if }\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{t}\end{array}\right.\quad \unicode[STIX]{x1D6FF}(r)=\left\{\begin{array}{@{}ll@{}}3\unicode[STIX]{x1D709}\quad & \text{if }\unicode[STIX]{x1D709}\leqslant \unicode[STIX]{x1D709}_{t}\\ h(r,t)/2\quad & \text{if }\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{t},\end{array}\right. & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D709}_{t}=h(T_{i}-T_{0})/6$ .

3.1.1 Thermal structure for an isoviscous flow, effect of $Pe$

The current cools by conduction and thermal boundary layers form at the contact with the surrounding medium. These boundary layers first connect at the tip of the flow, where the small thickness induces an important cooling (figure 3). A region of cold fluid forms at the front and slowly grows with time.

The flow is composed of a hot core, that contains a significant amount of heat and is lately referred to as the flow thermal anomaly, and a cold front. The radius $R_{c}(t)$ of the thermal anomaly is defined as the radius where $\unicode[STIX]{x1D6E9}_{b}=0.01$ . As the current thickens with time, a balance between advection and diffusion of heat is never reached in the current and the extent of the thermal anomaly grows with time. However, it spreads slower than the current itself and the extent of the cold fluid region at the tip, i.e. the region defined by $(R-R_{c})(t)$ , grows. For instance, for $Pe=100$ , while the region of cold fluid extends over approximately $10\,\%$ of the current at $t=0.5$ , it extends over approximately $20\,\%$ at $t=10$ (figure 3). The smaller the number $Pe$ , the more important the conductive cooling and the larger the cold region is (figure 4).

Figure 3.

Snapshots of the flow thermal structure $\unicode[STIX]{x1D703}(r,z,t)$ at different times indicated on the plot. Dashed lines represent the thermal boundary layers. Solid grey lines are isotherms for $\unicode[STIX]{x1D703}=0.2$ , $0.4$ , $0.6$ and $0.8$ . Here, the temperature field is decoupled from the dynamics, i.e. $\unicode[STIX]{x1D708}=1$ and $Pe=100$ .

Figure 4.

Snapshots of the flow thermal structure $\unicode[STIX]{x1D703}(r,z,t)$ for different set ( $\unicode[STIX]{x1D708}$ , $Pe$ ) with $\unicode[STIX]{x1D708}=1$ , $0.1$ , $0.01$ and $0.001$ and $Pe=1$ , $10$ , $100$ and $1000$ at $t=10$ . While $Pe$ controls the thermal structure of the flow, it has only a small influence on the flow aspect ratio which is controlled by $\unicode[STIX]{x1D708}$ .

Figure 5.

(a) Thickness normalized by the thickness at the centre $h(r,t)/h_{0}(t)$ versus radial axis normalized by the current radius $r/R(t)$ at different times indicated on the plot for $Pe=1.0$ and $\unicode[STIX]{x1D708}=1.0$ . The solid lines represent the thickness profiles and the dashed lines represent the thermal boundary layers. The predicted morphology (2.50) (dotted line) is also plotted for comparison. (b) Same plot but for $\unicode[STIX]{x1D708}=10^{-3}$ .

3.1.2 Thickness and temperature profile, effect of $\unicode[STIX]{x1D708}$

When accounting for the temperature dependence of the viscosity, the region of cold fluid at the tip is marked by a higher viscosity which enhances flow thickening at the expense of spreading. The larger the viscosity contrast, the larger the aspect ratio $h_{0}/R$ (figure 4). For instance, for the same value of $Pe=1$ , while the aspect ratio is $0.7$ for $\unicode[STIX]{x1D708}=1$ at $t=10$ , it is $4.2$ at the same time for $\unicode[STIX]{x1D708}=10^{-3}$ (figure 4). Nevertheless, the shape of the flow remains essentially self-similar, i.e. well described by (2.50) and cannot be differentiated from the shape of an isoviscous current if the thickness and the radial coordinates are rescaled by the thickness at the centre $h_{0}(t)$ and radius $R(t)$ (figure 5).

While the flow thermal structure is similar to the isoviscous case (figure 4), the important thickening induced by the viscosity increase tends to limit heat loss to the surrounding and to increase the size of the thermal anomaly at a given time. For instance, for $Pe=1$ at $t=10$ , while the thermal anomaly extends over approximately $30\,\%$ of the flow for $\unicode[STIX]{x1D708}=1$ , it extends over more than $50\,\%$ for $\unicode[STIX]{x1D708}=10^{-3}$ (figure 4).

As expected, a larger Péclet number leads to a larger thermal anomaly (figure 4). However, although different Péclet numbers cause very different thermal structures, the influence of the Péclet number on the flow morphology is small, much smaller than the effect of the viscosity contrast $\unicode[STIX]{x1D708}$ (figure 4). For instance, at $t=10$ for $\unicode[STIX]{x1D708}=10^{-3}$ , the thermal anomaly is still attached to the tip of the current for $Pe=1000$ whereas it makes approximately $50\,\%$ of the current for $Pe=1$ ; but, the thickness $h_{0}$ and the radius $R$ in both cases differ only by a few per cent (figure 4). This confirms that, in this regime, the spreading of the flow is not controlled by the mean temperature or average viscosity of the flow, but by the local dynamics at the flow front, as suggested by Lister et al. (Reference Lister, Peng and Neufeld2013).

3.2.1 Evolution of the thickness and the radius

In this bending dominated regime, the dynamics shows three different spreading phases. The thickness as well as the radius first evolve like an isoviscous flow, i.e. $h_{0}\propto t^{8/22}$ (2.52) and $R\propto t^{7/22}$ (2.53) (figure 6). In a second phase, thickening occurs at the expense of spreading and the dynamics deviates from the isoviscous case. Finally, the current returns to an isoviscous-like dynamics but offset from the hot isoviscous scaling law by a factor that depends on the viscosity contrast (figure 6).

Figure 6.

(a) Dimensionless thickness at the centre $h_{0}$ versus dimensionless time $t$ for different sets $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted lines: scaling laws $h_{0}=0.65h_{f}^{-1/11}\unicode[STIX]{x1D708}^{-2/11}t^{8/22}$ for $\unicode[STIX]{x1D708}=1.0$ and $0.01$ . (b) Dimensionless radius $R$ versus dimensionless time $t$ for the same sets $(\unicode[STIX]{x1D708},Pe)$ . Dotted lines: scaling laws $R=2.13h_{f}^{1/22}\unicode[STIX]{x1D708}^{1/11}t^{7/22}$ for $\unicode[STIX]{x1D708}=1.0$ and $0.01$ . (c) Dimensionless effective viscosity versus dimensionless time $t$ for different $Pe$ and $\unicode[STIX]{x1D708}=0.01$ . Solid lines: dimensionless effective viscosity $\unicode[STIX]{x1D702}_{e}$ defined by (3.3). Dotted lines: dimensionless average flow viscosity defined by $\overline{\unicode[STIX]{x1D702}_{a}(t)}=(1/V(t))\int _{0}^{R(t)}\int _{0}^{h(r,t)}r\unicode[STIX]{x1D702}(\unicode[STIX]{x1D703})\,\text{d}r\,\text{d}z$ where $V(t)$ is the current volume. Dashed lines: dimensionless average front viscosity $\unicode[STIX]{x1D702}_{f}$ defined by (3.4).

These different propagation phases thus reflect variations in the spreading rate of the current. As described in § 2.7, the propagation speed of the front is governed by local conditions in the peeling region and is given in dimensional form by (2.51); it decreases when the viscosity $\unicode[STIX]{x1D702}$ increases.

Inserting the cold viscosity $\unicode[STIX]{x1D702}_{c}$ in place of $\unicode[STIX]{x1D702}$ into (2.51), we find that the dimensionless thickness and radius for a cold isoviscous current evolve as

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle h_{0}=0.65\unicode[STIX]{x1D708}^{-2/11}h_{f}^{-1/11}t^{8/22}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle R=2.13\unicode[STIX]{x1D708}^{1/11}h_{f}^{1/22}t^{7/22}. & \displaystyle\end{eqnarray}$$

3.2.2 Flow effective viscosity

The effective viscosity of the current, i.e. the dimensionless viscosity $\unicode[STIX]{x1D702}_{e}$ which controls the current spreading rate in (2.51), is obtained by substituting $\unicode[STIX]{x1D708}$ by $1/\unicode[STIX]{x1D702}_{e}(t)$ in (3.1) and reads

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{e}(t)=\left(\frac{h_{0}(t)t^{-8/22}}{0.65h_{f}^{-1/11}}\right)^{11/2}, & & \displaystyle\end{eqnarray}$$

where $h_{0}(t)$ is given by our simulations.

As suggested by the results of § 3.2.1, the effective viscosity is first low (figure 6 c). It rapidly increases in the second phase of propagation and finally tends to the cold viscosity $1/\unicode[STIX]{x1D708}$ in the third phase. The effective viscosity is however much larger than the average flow viscosity (figure 6 c).

We calculate the average viscosity $\unicode[STIX]{x1D702}_{f}(t)$ over a front region of size $L$ in between $R(t)-L$ and $R(t)$

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{f}=\frac{1}{V_{f}}\int _{R-L}^{R}\int _{0}^{h}r\unicode[STIX]{x1D702}(\unicode[STIX]{x1D703})\,\text{d}r\,\text{d}z, & & \displaystyle\end{eqnarray}$$

where $V_{f}(t)$ is the volume of this region. The numerical evaluation of $\unicode[STIX]{x1D702}_{f}(t)$ for $L=4L_{p}(t)$ , where $L_{p}$ is the peeling length scale defined by Lister et al. (Reference Lister, Peng and Neufeld2013), given by

(3.5) $$\begin{eqnarray}\displaystyle L_{p}=1.08h_{f}^{13/22}\unicode[STIX]{x1D702}_{e}^{-2/11}t^{3/22}, & & \displaystyle\end{eqnarray}$$

fits relatively well the behaviour of the effective viscosity $\unicode[STIX]{x1D702}_{e}$ (figure 6 c).

Therefore, the effective viscosity, and thus the spreading rate, are controlled by the average viscosity of a small region at the front of size $L=O(L_{p})$ . Note that in the first flow phase, the effective viscosity $\unicode[STIX]{x1D702}_{e}$ , and thus the average viscosity of the peeling region, is larger than $1$ as the current has to accommodate the cold fluid initially present in the prewetting film. The effective viscosity is initially equal to 1 in the case of an initially hot prewetting film (see appendix C).

Average local thermal conditions in the peeling region thus control the spreading rate in the bending regime.

3.2.3 Evolution of the thermal condition at the tip

The thermal anomaly is first advected at the same velocity as the current itself, i.e. $R(t)=R_{c}(t)$ (figure 7 a). After a time that depends on $Pe$ and $\unicode[STIX]{x1D708}$ , the flow front leaves the thermal anomaly behind and $R(t)-R_{c}(t)$ increases with time (figure 7).

Figure 7.

(a) Extent of the cold fluid region $R(t)-R_{c}(t)$ versus dimensionless time for different combinations ( $\unicode[STIX]{x1D708}$ , $Pe$ ) indicated on the plot. (b) Same plot but where we rescale the extent of the cold fluid region by $Pe^{-1/3}\unicode[STIX]{x1D708}^{7/33}$ . Dotted line: scaling law $(R(t)-R_{c}(t))Pe^{1/3}\unicode[STIX]{x1D708}^{-7/33}=2.1h_{f}^{7/66}t^{9/22}$ .

In the bending regime, the interior pressure is constant and the thickness profile $h(r,t)$ is given by (2.50) (figure 5). The radius of the thermal anomaly $R_{c}(t)$ is theoretically taken as corresponding to the radius in the flow where heat advection locally balances heat loss, i.e.

(3.6) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}(\unicode[STIX]{x1D6E9}_{b}h)\sim Pe^{-1}\frac{\unicode[STIX]{x1D6E9}_{b}}{h}. & & \displaystyle\end{eqnarray}$$

Using the thickness profile (2.50), (3.6) becomes

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}^{2}\left(1+\frac{R_{c}}{R}\right)^{2}\left(\unicode[STIX]{x1D6E9}_{b}\frac{\text{d}h_{0}}{\text{d}t}+h_{0}\frac{\text{d}\unicode[STIX]{x1D6E9}_{b}}{\text{d}t}\right)+\frac{4h_{0}R_{c}^{2}\unicode[STIX]{x1D6E9}_{b}}{R^{3}}\frac{\text{d}R}{\text{d}t}\unicode[STIX]{x1D6FC}\left(1+\frac{R_{c}}{R}\right)\sim \frac{Pe^{-1}\unicode[STIX]{x1D6E9}_{b}}{\displaystyle \unicode[STIX]{x1D6FC}^{2}\left(1+\frac{R_{c}}{R}\right)^{2}h_{0}}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}(t)=(R(t)-R_{c}(t))/R(t)$ is the normalized region beyond $r=R_{c}(t)$ . In the limit $\unicode[STIX]{x1D6FC}\ll 1$ , i.e. $R_{c}/R\sim 1$ , the time derivative is locally dominated by its advective part ( $\propto \unicode[STIX]{x1D6FC}$ ) and we finally get

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}^{3}\sim \frac{Pe^{-1}}{h_{0}^{2}(t)}\frac{R}{\displaystyle \frac{\unicode[STIX]{x2202}R}{\unicode[STIX]{x2202}t}}. & & \displaystyle\end{eqnarray}$$

Substituting $h_{0}(t)$ and $R(t)$ by their respective scaling laws (3.1) and (3.2), the size evolution of the normalized cold front region $\unicode[STIX]{x1D6FC}$ reads

(3.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(t)\sim Pe^{-1/3}\unicode[STIX]{x1D708}^{4/33}h_{f}^{2/33}t^{1/11}, & & \displaystyle\end{eqnarray}$$

which is equivalent to

(3.10) $$\begin{eqnarray}\displaystyle R(t)-R_{c}(t)\sim Pe^{-1/3}\unicode[STIX]{x1D708}^{7/33}h_{f}^{7/66}t^{9/22}. & & \displaystyle\end{eqnarray}$$

The predicted scaling law for the evolution of the cold fluid region (3.10) closely fits the numerical simulations for $\unicode[STIX]{x1D708}<1$ and for different Péclet numbers (figure 7). In particular, this scaling perfectly fits the simulations with a numerical prefactor equal to $2.1$ .

3.2.4 Summary of the bending regime dynamics

The spreading rate of the current is controlled by the average viscosity in the peeling region which depends on the average local thermal conditions at the tip.

At the initiation of the flow, heat advection is larger than conductive heat losses in the peeling region, the thermal anomaly reaches the flow front and the spreading rate is similar to the one of an isoviscous hot current.

Once heat loss in the peeling region becomes larger than heat advection, the current tip leaves the thermal anomaly behind. The average temperature of the peeling region decreases and the effective viscosity rapidly increases. Setting $(R-R_{c})(t)$ (3.10) equal to the peeling length scale $L_{p}(t)$ (3.5), taking $\unicode[STIX]{x1D708}=1$ and inverting for time, we found that the time $t_{b2}$ to reach the second thermal phase reads

(3.11) $$\begin{eqnarray}\displaystyle t_{b2}=7\times 10^{-3}Pe^{11/9}h_{f}^{16/9}, & & \displaystyle\end{eqnarray}$$

where the numerical prefactor is chosen from the numerical solutions. When rescaling the time of the simulations by $t_{b2}$ , the different simulations enter the second phase simultaneously (figure 8 a).

Figure 8.

(a) Dimensionless thickness at the centre $h_{0}/h_{0}(t_{b2})$ versus dimensionless time $t/t_{b2}$ , where $t_{b2}$ is the time to enter the second thermal phase of the bending regime (3.11), for different sets $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted line: scaling laws $h_{0}(t)/h_{0}(t_{b2})=(t/t_{b2})^{8/22}$ . (b) Dimensionless thickness at the centre $h_{0}/h_{0}(t_{b3})$ versus dimensionless time $t/t_{b3}$ , where $t_{b3}$ is the time to enter the third thermal phase of the bending regime (3.12), for different sets $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted line: scaling laws $h_{0}(t)/h_{0}(t_{b3})=(t/t_{b3})^{8/22}$ .

Finally, when the peeling region becomes entirely cold, i.e. $R(t)-R_{c}(t)\gg L_{p}(t)$ , the flow behaves as an isoviscous cold current. Repeating the same exercise as before while keeping the viscosity contrast $\unicode[STIX]{x1D708}$ , we found that the time $t_{b3}$ for the flow to enter this third phase reads

(3.12) $$\begin{eqnarray}\displaystyle t_{b3}=14Pe^{11/9}h_{f}^{16/9}\unicode[STIX]{x1D708}^{-1/9}, & & \displaystyle\end{eqnarray}$$

where the numerical prefactor is chosen from the numerical solutions. When rescaling the time of the simulations by $t_{b3}$ , the different simulations enter the third phase simultaneously (figure 8 b).

4.1.1 Thermal structure for an isoviscous flow, effect of $Pe$

As in the bending regime, the bulk of the fluid first expands at the injection temperature and $R_{c}(t)=R(t)$ . As the bottom and top cool by conduction, thermal boundary layers form at the contact with the surrounding medium and connect at the tip of the current. However, in the gravity current regime, for a constant viscosity, the thickness of the current tends to a constant. Therefore, when conduction in the surrounding medium balances the input of heat at the centre and the flow front has already left the thermal anomaly behind, its extent approaches a steady state (figure 9).

Figure 9.

Snapshots of the flow thermal structure $\unicode[STIX]{x1D703}(r,z,t)$ at different times indicated on the plot. Dashed lines: thermal boundary layers. Here, $\unicode[STIX]{x1D708}=1$ and $Pe=100$ .

The radius of the steady-state thermal anomaly $R_{c}$ also largely depends on $Pe$ in this regime: the larger the number $Pe$ , the larger the radius $R_{c}$ (figure 10).

Figure 10.

Snapshots of the flow thermal structure $\unicode[STIX]{x1D703}(r,z,t)$ for different sets ( $\unicode[STIX]{x1D708}$ , $Pe$ ) with $\unicode[STIX]{x1D708}=1$ , $0.1$ , $0.01$ and $0.001$ and $Pe=1$ , $10$ , $100$ and $1000$ at $t=200$ .

4.1.2 Thickness and temperature profile, effect of $\unicode[STIX]{x1D708}$

For a current with a viscosity that depends on temperature, as soon as the cooling becomes more effective and the thermal anomaly detaches from the current tip, the spreading slows down and flow thickening is enhanced (figure 10). For instance, for $Pe=1$ , while the aspect ratio $h_{0}/R$ is approximately $0.12$ for $\unicode[STIX]{x1D708}=1$ at $t=200$ , it is ${\sim}1$ for $\unicode[STIX]{x1D708}=10^{-3}$ (figure 10). The shape of the current is not self-similar and the front steepens when the viscosity increases in comparison to the isoviscous case as noted by Bercovici (Reference Bercovici1994). However, when the current becomes much larger than the thermal anomaly, the current side slumps to become less steep (figure 10) and recovers a shape similar to an isoviscous flow with a cold viscosity.

The thermal structure is similar to the isoviscous case ( $\unicode[STIX]{x1D708}=1$ ). In particular, after a time that depends on $Pe$ , the thermal anomaly approaches a steady-state profile (figure 10). As in the bending regime, the thickening at the centre limits heat loss to the surrounding for large values of the viscosity contrast $\unicode[STIX]{x1D708}$ . Therefore, the extent of the thermal anomaly in the steady state is slightly larger for a larger viscosity contrast. For instance, for $Pe=10$ at $t=200$ , while the thermal anomaly extends over less than $2$ for $\unicode[STIX]{x1D708}=1$ , it reaches ${\sim}3$ for $\unicode[STIX]{x1D708}=10^{-3}$ .

The flow morphology is much more sensitive to $Pe$ in the gravity current regime than in the bending regime and different $Pe$ lead to different current morphologies for a given $\unicode[STIX]{x1D708}$ (figure 10). For instance, for $\unicode[STIX]{x1D708}=10^{-3}$ at $t=200$ , the thermal anomaly is still attached to the tip for $Pe=10^{3}$ and the aspect ratio of the flow $h_{0}/R$ is close to $0.15$ . In contrast, for $Pe=1$ , the thermal anomaly radius is less than $30\,\%$ of the current radius and the dimensionless aspect ratio of the flow is much larger $h_{0}/R=1.15$ (figure 10).

4.2.1 Evolution of the thickness and radius

As in the bending regime, the dynamics in the gravity current regime shows three different spreading phases. The thickness as well as the radius first follow the isoviscous scaling laws for a given hot viscosity $\unicode[STIX]{x1D702}_{h}$ , i.e. $h_{0}$ approaches a constant and $R\propto t^{1/2}$ (figure 11). In a second phase, the thickness rapidly increases and the spreading slows down. Finally, the current returns to an isoviscous-like dynamics but offset from the hot isoviscous scaling law by a factor that depends on the viscosity contrast (figure 11). In particular, replacing $\unicode[STIX]{x1D702}$ by $\unicode[STIX]{x1D702}_{c}$ instead of $\unicode[STIX]{x1D702}_{h}$ in (2.54), we find that the dimensionless radius $R(t)$ in the third flow phase matches the one of a cold viscosity current (4.1) (figure 11)

(4.1) $$\begin{eqnarray}\displaystyle R(t)=1.10\unicode[STIX]{x1D708}^{1/8}t^{1/2}. & & \displaystyle\end{eqnarray}$$

As in the bending regime, these spreading rate variations reflect variations in the flow effective viscosity, from the hot to the cold viscosity. However, in that case, the effective viscosity is the average flow viscosity. Using the predicted scaling law for the radius $R(t)$ (4.1) as a function of $\unicode[STIX]{x1D708}_{e}$ to derive the effective viscosity $1/\unicode[STIX]{x1D708}_{e}$ , we show that the current average viscosity matches the effective flow viscosity almost perfectly (figure 11 c).

Figure 11.

(a) Dimensionless thickness at the centre $h_{0}$ versus dimensionless time $t$ for different sets $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted lines represent the scaling laws $h_{0}=2.1\unicode[STIX]{x1D708}^{-1/4}$ for $\unicode[STIX]{x1D708}=1.0$ and $10^{-2}$ . (b) Dimensionless radius $R$ versus dimensionless time $t$ for the same sets $(\unicode[STIX]{x1D708},Pe)$ . Dotted lines represent the scaling laws $R=1.1\unicode[STIX]{x1D708}^{1/8}t^{1/2}$ for $\unicode[STIX]{x1D708}=1.0$ and $10^{-2}$ . (c) Dimensionless effective viscosity versus dimensionless time $t$ for different $Pe$ and $\unicode[STIX]{x1D708}=0.01$ . Solid lines: effective viscosity $\unicode[STIX]{x1D702}_{e}(t)=(R(t)/1.1t^{1/2})^{-8}$ . Dashed lines: average flow viscosity defined by $\overline{\unicode[STIX]{x1D702}_{a}(t)}=(1/V(t))\int _{0}^{R(t)}\int _{0}^{h(r,t)}r\unicode[STIX]{x1D702}(\unicode[STIX]{x1D703})\,\text{d}r\,\text{d}z$ where $V(t)$ is the current volume. Note that the irregularities at early times, which are enhanced by the log–log representation, come from a slight decrease of the effective viscosity when the current begins to thickens.

4.2.2 Characterization of the thermal anomaly

The thermal anomaly is first advected at the same velocity as the current itself, i.e. $R_{c}(t)/R(t)=1$ (figure 12 a). After a time that depends on $Pe$ and $\unicode[STIX]{x1D708}$ , the flow front leaves the thermal anomaly behind and the anomaly extent slowly approaches a steady-state profile as the viscosity becomes closer to the cold viscosity (figures 9, 10 and 12).

Figure 12.

(a) Normalized thermal anomaly radius $R_{c}(t)/R(t)$ versus dimensionless time for different combinations ( $\unicode[STIX]{x1D708}$ , $Pe$ ) indicated on the plot. (b) Same plot but where we rescale the normalized thermal anomaly radius $R_{c}(t)/R(t)$ by $Pe^{1/2}\unicode[STIX]{x1D708}^{-1/4}$ . Dotted line: $Pe^{-1/2}\unicode[STIX]{x1D708}^{1/4}R_{c}(t)/R(t)=0.8t^{-1/2}$ .

At the steady-state radius $R_{c}$ of the thermal anomaly, a balance between heat advection and diffusion in the surrounding medium gives

(4.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}_{b}\frac{U_{0}}{R_{c}}\sim \frac{Pe^{-1}}{h_{0}^{2}}\unicode[STIX]{x1D6E9}_{b}, & & \displaystyle\end{eqnarray}$$

where $U_{0}$ is a mean dimensionless velocity of advection. For a gravity current, in contrast to the bending regime, the thickness $h_{0}$ approaches a constant. Taking $U_{0}$ as a horizontal redistribution of the dimensionless injection rate at $r=R_{c}$ , i.e. $U_{0}=1/(R_{c}h_{0})$ , and using the scaling for the thickness of a cold isoviscous current $h_{0}\sim \unicode[STIX]{x1D708}^{-1/4}$ , we obtain

(4.3) $$\begin{eqnarray}\displaystyle R_{c}\sim Pe^{1/2}\unicode[STIX]{x1D708}^{-1/8} & & \displaystyle\end{eqnarray}$$

and hence

(4.4) $$\begin{eqnarray}\displaystyle \frac{R_{c}}{R(t)}\sim Pe^{1/2}\unicode[STIX]{x1D708}^{-1/4}t^{-1/2}, & & \displaystyle\end{eqnarray}$$

where we have used (4.1).

This scaling law fits our simulation for a prefactor equal to $0.8$ (figure 12): when the thermal anomaly enters the steady state, its radius remains constant and the normalized radius $R_{c}(t)/R(t)$ evolves as the inverse of the current radius, i.e. as $t^{-1/2}$ (figure 12). Furthermore, both the dependence with $Pe$ and $\unicode[STIX]{x1D708}$ vanish when rescaling $R_{c}/R(t)$ by $Pe^{1/2}\unicode[STIX]{x1D708}^{-1/4}$ in the steady state (figure 12 b).

4.2.3 Summary of the gravity regime dynamics

In the gravity regime, the current spreading rate is controlled by the flow average viscosity. Therefore, it strongly depends on the extent of the thermal anomaly.

At flow initiation, the thermal anomaly is advected at the same velocity as the current itself and the current spreads with a hot viscosity $\unicode[STIX]{x1D702}_{h}$ . Once heat loss over the current balances the heat input at the centre, the current tip leaves the thermal anomaly behind, the average viscosity increases and the spreading rate decreases. The time $t_{g2}$ to enter this second phase scales with the time to cool the hot current ( $\unicode[STIX]{x1D708}=1$ ) by conduction and reads

(4.5) $$\begin{eqnarray}\displaystyle t_{g2}=0.65Pe, & & \displaystyle\end{eqnarray}$$

where the numerical prefactor is chosen from the numerical solutions. Indeed, when rescaling the time of the simulations by $t_{g2}$ , the different simulations enter the second phase simultaneously (figure 13 a).

Figure 13.

(a) Dimensionless radius at the centre $R/R(t_{g2})$ versus dimensionless time $t/t_{g2}$ , where $t_{g2}$ (4.5) is the time to enter the second thermal phase of the gravity regime, for different sets of $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted line: scaling law $R/R(t_{g2})=(t/t_{g2})^{1/2}$ . (b) Dimensionless radius $R/R(t_{g3})$ versus dimensionless time $t/t_{g3}$ , where $t_{g3}$ (4.6) is the time to enter the third thermal phase of the gravity regime, for different sets $(\unicode[STIX]{x1D708},Pe)$ indicated on the plot. Dotted line: scaling law $R/R(t_{g3})=(t/t_{g3})^{1/2}$ .

Finally, when the thermal anomaly becomes small compared to the current, i.e. $R_{c}/R\ll 1$ , the average flow temperature is close to zero and the current behaves as an isoviscous cold current. The time $t_{g3}$ to enter this third phase scales as the time to cool an isoviscous cold gravity current and reads

(4.6) $$\begin{eqnarray}\displaystyle t_{g3}=10Pe\unicode[STIX]{x1D708}^{-1/2}, & & \displaystyle\end{eqnarray}$$

where the numerical prefactor is chosen from the numerical solutions. Indeed, when rescaling the time of the simulations by $t_{g3}$ , the different simulations enter the third phase simultaneously (figure 13 b). The time scales $t_{g2}$ (respectively $t_{g3}$ ) can be derived by inverting (4.4) for $t$ for a fixed value of $R_{c}/R$ and setting $\unicode[STIX]{x1D708}=1$ (respectively keeping $\unicode[STIX]{x1D708}$ ) in a similar way that we derive $t_{b2}$ (respectively $t_{b3}$ ) in § 3.2.4.