Non-isothermal thin-film flow of a viscoplastic material over topography: critical Bingham number for a partial slump

This paper considers the non-isothermal flow of a viscoplastic fluid on a horizontal or an inclined surface with a flat, a step-up and a step-down topography. A particular application of interest is the spread of a fixed mass—a block—of material under its own weight. The rheology of the fluid is described by the Bingham model which includes the effect of yield stress, i.e. a threshold stress which must be exceeded before flow can occur. Both the plastic viscosity and the yield stress are modelled with temperature-dependent parameters. The flow is described by a reduced model with a thin-film equation for the height of the block and a depth-averaged energy conservation equation for the heat transfer. Results show that for large values of the yield stress, only the outer fraction of the fluid spreads outward, the inner fraction remaining unyielded, hence the block only partially slumps. Conversely, for small values of the yield stress, the entire block of fluid becomes yielded and therefore slumps. We present an analysis which predicts the critical value of the yield stress for which partial slump occurs and how it depends on temperature.


Introduction
As the name suggests, a viscoplastic fluid can either behave as a viscous fluid, which deforms continuously when subject to a sufficiently strong mechanical load, or as a solid, i.e. the material does not deform. As an illustrative example, a volume of viscous fluid deposited on a flat surface will spread under its own weight until the gravity and capillary forces are at equilibrium. A viscoplastic fluid, on the other hand, will only start spreading if the gravity force is strong enough to overcome the fluid cohesive force, this explains why a lump of toothpaste rests on a surface without spreading. The flow/no flow boundary for a viscoplastic fluid is dictated by the yield stress, a rheological property of the fluid. Flow can only occur when the stress within the fluid exceeds the yield stress. Viscoplastic fluids appear in a wide range of contexts including food products (e.g. ketchup, chocolate, peanut butter), healthcare products (e.g. toothpaste, creams), the petroleum industry, cosmetics, geophysical flows (lava is known to have a viscoplastic behaviour [3,19,22,34]), and many others.
Three-dimensional and two-dimensional isothermal thin-film flows, such as gravity-driven flow over an inclined plane, spreading flows and spin coating, are some of the fundamental problems in fluid mechanics. Thin-film flow of a viscoplastic material, on a surface either flat or with topographical features, has been extensively studied in the past, both theoretically and numerically, see [4][5][6]11,20,21,25,26,33,39,43] and references therein. Experimentally, these kind of set-ups have been used to infer rheological properties of the material, such as the yield stress. For an excellent review of experimental data involving yield stress fluids see [14]. For example, in the case of gravity-driven flow over an inclined plane, it has been shown that the critical thickness associated with stoppage for a wide uniform layer of a yield stress fluid is directly proportional to the yield stress, see [15,16]. Similarly, spreading flows can be used to predict the yield stress, whereby if the initial volume poured is enough such that the critical thickness is exceeded, then the material flows, see [17]. For the spin coating set-up, Tabuteau et al. [42], showed that a sample of viscoplastic material starts to flow beyond a critical angular velocity.
We should note that all the experimental, theoretical and numerical results stated above were done for isothermal conditions. To our knowledge, equivalent results for non-isothermal flows are lacking in the literature. How does temperature affect the critical thickness of a slump of viscoplastic material? Is there a temperature-dependent critical yield stress for which a partial slump exists? These are some of the questions we try to answer in this work.
When considering non-isothermal dynamics, one has to take into account that as the hot material spreads (or flows), it cools down by radiation and convection to the surroundings and by conduction to the ground, and that many rheological parameters are temperature-dependent. The effects of temperature-dependent material properties on the spreading of gravity currents have been considered by numerous authors. For example, the effects of a temperature-dependent viscosity were included in the lubrication approximation by Hutter [24], Oron et al. [35], Wilson and Duffy [45] and Sansom et al. [38] to investigate the draining of a non-isothermal thin liquid film on a surface. The latter approach was used by Boujo and Sellier [10] in the context of the control of non-isothermal gravity currents. Vasilyev et al. [44] shed light on the effects of viscous dissipation upon gravity currents while Osiptsov [36] derived a family of hydrodynamic models of a cooling lava flow over a conical surface. Lyman et al. [31] studied experimentally the effects of surface solidification and yield strength on a fixed volume of fluid released in a horizontal channel. Zadražil et al. [46] considered the spreading and imbibition dynamics of hot two-dimensional droplets on a cooler porous medium. Garel et al. [18] developed a thermal model for the cooling of an axisymmetric isoviscous gravity current. The model includes surface cooling by radiation and convection, and basal cooling by conduction to a substrate. Recently, Pascal et al. [37] studied the stability of a power-law fluid film flowing down a heated incline. For an extensive review on non-isothermal rheology of viscoplastic materials see [12].
A non-isothermal thin-layer theory for a viscoplastic material was first developed by Balmforth and Craster in [1]. There, they study axysimmetric viscoplastic domes, with constant extrusion velocity, described by the Herschel-Bulkley constitutive law. In order to derive an energy equation consistent with a lubrication approximation the authors keep only zeroth and second-order terms. The first-order approximation is assumed to be zero and the leading order one does not depend on the vertical direction. The resultant cooling law is a balance between a diffusion term which is O(Pe −1 ), here Pe is the Péclet number (ratio between convection and diffusion), and a linear cooling law at the surface which is O(h −1 ), with h being the height of the material. Because of the asymptotic nature of the model, the validity is restricted to Pe ∼ O(1), hence the material always cooled. This problem was fixed by the same authors in [2] where they consider the full energy equation, i.e. no averaging across the vertical direction, and solve the problem numerically; they called this model the "Full thin-layer model". Their results show that for large values of Pe, cooling is concentrated at the upper surface of the fluid and a chilled boundary layer descends into the interior of the fluid as it expands. They found that temperature can be described surprisingly well by a polynomial in the vertical direction, corroborating the experimental results of [40] and the modelling work of [7,8,18,30]. A second degree polynomial was chosen following the work in [8], where it was stated that even though higher order contributions to the vertical temperature profile may be important, the parabolic contribution must always be the dominating term, and the higher-order contributions typically decay away relatively rapidly because of their large temperature gradient. The authors then develop a skin-theory model for the boundary layer and compare results using the full thin-layer model, the skin-theory model and the vertically isothermal model of [1]. In the absence of yield stress and with a temperature-dependent viscosity, the models compare well for Pe ≤ 1. On the other hand, if Pe is large, the vertically isothermal model is inadequate and the results deviate completely from the other two models, which compare really well. Interestingly, all three models describe the spreading dynamics of the domes relatively well if the yield stress is different from zero. It is shown that thermal effects are damped by the existence of a pseudo-plug. In a more recent work, Bernabeu et al. [9] combined the approaches of Bercovici et al. [8], and Balmforth et al. [1,2], and developed a reduced model consisting of a thin-film equation for the viscoplastic material and an averaged heat equation for the temperature. They assume a third degree polynomial description for the temperature in the vertical direction. They tested their approach against the assumption of a parabolic profile in the vertical direction and they found that the difference in results is negligible. The investigation of a temperature-dependent critical yield stress for which a partial slump exists was not considered in the references above.
In this paper, we build on the work of Bernabeu et al. [9], we develop and implement a mathematical model describing the flow of a non-isothermal viscoplastic fluid either spreading under its own weight on a horizontal surface or draining down an inclined plane. Of particular interest here is the effect of heat transfer on the spreading and rest state of a viscoplastic slump.
The rest of the paper is organised as follows. In Sect. 2 we present the full model and its nondimensionalisation. In Sect. 3 we derive the reduced non-isothermal viscoplastic model, following the work in [9]. In Sect. 4 we discuss the numerical method used to solve the evolution equations for the height of the material and temperature profile. In Sect. 5 we present the results of our study. Firstly, we consider a horizontal flat topography with no extrusion velocity in Sect. 5.1, and under these conditions we derive a semi-analytical expression for the critical Bingham number for which a partial slump occurs in Sect. 5.1.1. Secondly, in Sect. 5.2 we extend previous results for the cases of step-up and step-down topographies. Thirdly, in Sect. 5.3 we consider the case of a nonzero extrusion velocity on a horizontal flat surface and present spreading dynamics of the block and decay rate of the temperature as a function of extrusion velocity. Finally, some concluding remarks are provided in Sect. 6.

Problem description
We will consider the flow of a block of material, described by a viscoplastic constitutive relation, under the effects of gravity over a topography described by a functional relation. The main goal of our work is to study thermal effects on the flow, hence a temperature-dependent rheology is to be considered. Assuming a Boussinesq approximation, the conservation equations describing the flow are: where D/Dt is the material derivative, u = (u, v, w) is the velocity field, p is the pressure, τ is the stress tensor, g is the gravity, is the temperature, κ is the thermal diffusivity, ρ and c p are the density and specific heat of the fluid. Note that we have assumed that viscous dissipation effects are negligible in Eq. (3). The domain is specified in the following way: where ⊂ R 2 is a region in the plane, h(t, x, y) is the height of the free surface, f (x, y) is the functional representation of the topography, as shown in Fig. 1. Boundary conditions for Eqs. (1) and (2) are no-slip at the bottom and stress free at the free surface. For Eq. (3) we will consider no heat transfer at the bottom and convection fluxes at the free surface.
In Eq. (1) the deviatoric stress τ is defined as: where μ( ) and τ y ( ) are the temperature-dependent plastic viscosity and yield stress, respectively. The second invariants of the stress and rate of strain tensors are defined as: Assuming that the characteristic length of the height of the material is much smaller than the characteristic length of the terrain, we consider the dimensionless aspect ratio where H is the vertical scale and L the horizontal scale. We non-dimensionalise our spatial variables using with velocity, pressure, temperature and temporal scales where b is the characteristic temperature which is set to the initial temperature of the block or the temperature of the extruded material, and μ p is the viscosity of the material at temperature b . Substituting into Eqs. (1) to (3) and dropping the hats for the non-dimensional variables, we get: Here, and Pe = ρc p U L κ are the Reynolds and Péclet number, respectively. We will assume Re ∼ O(1) and ε > Pe −1 , hence Pe 1. In particular, we will definePe = ε 2 Pe.
The non-dimensional constitutive relation is given by: where the temperature-dependent Bingham number is defined as: where τ y is the yield-stress of the material at temperature b . Two points to note in expression (12), as we are interested in temperature-dependent rheology, in this case the yield stress, f B ( ) is a non-dimensional function of temperature to be determined. The second point is that we will assume B i ∼ O(1), this means that the ratio between yield stress and gravity effects is small, i.e. the material will always deform.

Reduced model
The main assumptions for the derivation of the reduced model are: • Small Reynolds number, i.e. creeping flow.
• Topographical features with long wavelength and small slope.
• Small aspect ratio between the vertical and horizontal directions.
• Second degree polynomial as an ansatz for temperature-dependence in the vertical direction.
As it is standard in lubrication theory, we consider the O(1) approximation of system (6)- (10): with boundary conditions We now derive the lubrication equations by averaging our variables in the vertical direction from z = f (x, y) to z = f (x, y) + h(t, x, y). Integrating Eqs. (13)-(16) we get the following evolution equation for the free surface: where the ∇-operator takes the gradient in the (x, y)-direction only, w s is the depth-averaged vertical velocity which in this case acts as a source term, and h c is the position of the yield surface of the material given by: Note that for simplicity we have made the assumption that the rheology depends only on the depth-averaged temperature: In order to get the reduced energy equation, we take the depth-average of Eq. (17): where the overbar notation means averaging in the z-direction, and u xy = (u, v).
In order to deal with the second term of the left hand side of (20 and at the free surface and at the bottom in which by doing this Eq. (20) becomes It is necessary to add an extra closure relationship if we want to avoid having to calculate the integral term ( ) in Eq. (24) in our numerical simulations. In [9], the authors assumed φ(t, x, y, z) to be a cubic polynomial in z and included the closure relationship φu xy = u xy (25) such that the equality in (25) holds. We assume here for Eq. (25) to be zero while keeping the quadratic polynomial assumption for φ. The reasoning is as follows: (i) In [9], the authors performed numerical simulations using both a cubic and a quadratic polynomial and found the difference between the results to be negligible; (ii) The extra added coefficient, a 3 (t, x, y) only affects the rate at which the material cools and is smaller in magnitude than a 2 (t, x, y). Therefore, Eq. (24) simplifies to For simplicity we will consider that topographical changes only occur in the x-direction of the flow, Fig. 1. This reduces our model to an unsteady one-dimensional model for the free surface position and the averaged temperature. To avoid any notational confusion we will drop the "av" subscript in the temperature value and from the rest of the paper (t, x) corresponds to the vertically averaged temperature of the fluid. The system of partial differential equations to solve is: where the yield-surface is given by: The rheological parameters depend on temperature in the following way, see [9]: where α μ and α B are non-dimensional constants that control the exponential decay for the temperature dependent rheological parameters. For Eq. (28), the averaged velocity is given by the following expression, The coefficient a 2 (t, x) in Eq. (28) is the pointwise solution of a 3 × 3 system of equations that depend on the conditions (21) to (23). We can eliminate one of the equations and reduce the problem to a 2 × 2 system. Therefore, in order to find a 2 (t, x), we need to solve: Note that the spatial domain is the position of the left advancing front and x fr (t) is the position of the right advancing front. We impose homogenous Dirichlet boundary conditions at the fronts for both h(t, x) and (t, x). Initial conditions are in which we use r (x) for both h init and init . We define the topography using the functional relation where m is the slope of the inclined defined as m = tan(θ ) (θ is the angle of inclination defined in Fig. 1), a is the amplitude (height) of the steps, b is a constant that controls the location of the step and c is related to the steepness of the step. Note that by using this definition of f (x) we have the flexibility to study inclined or horizontal planes with flat, step-up and step-down domains.

Numerical method
To avoid having to compute the solution of our system in a very large domain, we map our moving physical domain to a fixed computational domain and advance our left and right fronts through some evolution equations. This technique has been applied for viscoplastic dam breaks in [4,32]. We let where x fl (t) and x fr (t) are the left and right fronts, respectively. To simplify our notation we will define Substituting Eq. (36) into Eqs. (27) and (28) we get the system: where() means time differentiation. Now We have used an implicit-explicit first-order method for time integration. Expressions (33), (39) and (40) are treated explicitly. Second-order derivatives in space are treated implicitly. We use a second-order central difference scheme for Eq. (37). Equation (38) is hyperbolic. In order to get a stable scheme we use a first order upwind scheme. For all our simulations, we set x = 10 −2 and t = 10 −3 . In order to get evolution equations for x fl and x fr , we set h t (t, x fl ) = h t (t, x fr ) = 0 ∀t in Eq. (39). Clearly, when x = x fl we have ξ = 0 and we get the following evolution equation for x fl : Similarly, the evolution equation for x fr reads: To integrate Eqs. (41) and (42) we use an explicit Euler scheme with the same time step as before. We run our simulations until x b reaches a steady state by taking the difference between two successive time steps and checking that it is less than a set tolerance which for our work is 10 −7 .
In the following section we will compare results between isothermal and non-isothermal models.

Results
In this section we investigate the effects that a temperature-dependent rheology has on the spread rate, length of the block, shape of the free surface and critical yield stress. We will compare these results with the isothermal model. Several topographies will be considered. In Sects. 5.1 and 5.2, we consider w s (t, x) = 0. As discussed earlier, even though we are interested in advection dominated flows, in our lubrication approximation our modified Péclet numberPe could be several orders of magnitude smaller than the bulk Péclet number and diffusive effects play a role in the flow. Red and blue curves are the non-isothermal cases and the black broken curves are the isothermal case. We consider two values forPe, red curves correspond toPe = 1000 and blue curves correspond toPe = 100. First thing to note is that for the smaller Péclet number problem, the temperature has already decayed to zero for t = 1000. As we can see, even though the decay of the temperature is significant, the isothermal and non-isothermal height profiles do not differ much, especially for the viscoplastic case. This is confirmed in Fig. 3 where we show the length of the block as a function of time for the three Bingham numbers. The Newtonian block spreads as t n with n = 0.197 (isothermal) and n = 0.188 (non-isothermal cases). The decay rate for the isothermal case corresponds to the theoretical prediction of n = 1/5, see [28].
As it was shown by [32], when yield stress is present, there exists an arrested state for which the block stops deforming, confirmed in Fig. 2. One of the characteristics of viscoplastic materials is that for large values of Bingham numbers there exist a yield point x y such that for all |x| ≤ x y the material does not deform. For the case of B i = 1 we have |x y | = 0.72, see Fig. 2c. As we can see in Fig. 3 the final length of the block for the viscoplastic material decreases asPe decreases. Clearly, the temperature-dependent yield stress affects the deformation and spreading of the material and we would expect that the critical Bingham number, that is the minimum Bingham number for which a yield point exist, must be affected by the temperature. We will explore this dependence in the next section.

Critical Bingham number, B i,c , semi-analytical results
As we mentioned before, the main characteristic of the spreading of a viscoplastic material is that it attains an arrested state, i.e. h c → 0 as t → ∞. This arrested state, h ∞ (x), can be found when the second argument in Eq. (29) is exactly zero, i.e. when Note that by satisfying the balance in Eq. (43), the yield surface is equal to zero over the whole domain and we have no flow. When considering a horizontal flat topography, an analytical solution exists for the steady-state profile h ∞ (x), see for example [4,32] Here, x f,∞ is the position of the front of the slump and it is defined as 2. For a partial slump we get where x y,∞ is the yield point. Both yield point and front are given by the following expressions Expression (47) gives us the required condition to find the critical Bingham number. If B i,c = 2/3 then x y,∞ = 1/2.
The goal of this section is to investigate the effect that a temperature-dependent yield stress has on the value of the critical Bingham number. Clearly, this will depend on the functional description of f B ( ). As we have stated in Eq. (31), we will consider an exponential law for the yield stress. Now, we are in a position to investigate the decay rate of as t → ∞.
From Eq. (28) we can write the characteristic equations Solving system (33) for the case f = 0 we get .
Equations (49) and (50) can be integrated to get Note thatū = 0 for x = 1/2 for all t and because of symmetry is maximum there, hence In order to investigate the decay rate of Eq. (53), we solve system (37) and (38) for various values ofPe, α B and B i and track (t, 1/2). For each time series we calculate log(|| || ∞ ) and plot against t. We found that where C(Pe, B i ) is a constant that depends only onPe and B i . This suggests that the integrand in Eq. (53) is independent of t. If we consider then with We now assume that Eq. (53) decays as exp(−A(B i )t/Pe). In Fig. 4 we plot −Pe log(|| || ∞ )/A(B i ) against t for different values ofPe, B i and α B . The data collapses into the line y(t) = t, confirming our assumption. Because ∼ O (1) and in general α B ≤ 1, a good approximation to Eq. (31) is Backed by our numerical simulations we found that for t ∼ O(1) the block has mostly spread and the maximum of h does not move much, hence Substituting Eq. (59) into Eq. (58) we get an expression for the temperature-dependent Bingham number for values close to its critical value Using Eq. (60) in Eq. (47) we can get an expression for the critical Bingham number by solving In Fig. 5, we present results of the critical Bingham number, B i,c , as a function of α B andPe. In Fig. 5a we consider the full problem, i.e. we solve system (37) and (38) for fixedPe and α B until we reach steady state and use the bisection method to find B i such that h(1/2) = 1. We choose the bisection method because of its global convergence properties but the main drawback is its sublinear convergence rate. Even with a good initial interval, one needs to take about 8 iterations for each pair (α B ,Pe). This is a very slow process and also, the value of B i,c is dependent on the accuracy of the numerical method chosen. On the other hand, Fig. 5b shows the solution of Eq. (61) for a given pair (α B ,Pe). This can be done using a standard root finding algorithm (in our case we used fzero in MATLAB). In Fig. 5c we plot slices of the surfaces in Figs. 5a and b for α B = 0.5 and 1. As we can see, our approximation works well if α B /Pe is small. For larger values of the ratio, the analytical solution over-estimates the value for B i,c . A possible solution could be to modify Eq. (59) in order to capture the correct rate at which the centre of the block cools. One feature of the two approaches that is worth mentioning is that as the ratio α B /Pe gets larger, thermal effects strongly affect the value of the critical Bingham number. This makes sense asPe dictates the rate at which the temperature goes to zero and α B controls the non-isothermal effects in the rheology of the material. Therefore, as we can see in Eq. (61), for the case of small α μ , the critical Bingham number depends mainly on the ratio α B /Pe. In Fig. 5d, we show B i,c as a function of α B /Pe. Clearly, as α B /Pe → ∞ then B i,c ≈ K (α B /Pe) −n , with K a constant of proportionality. The full black curve corresponds to for α B /Pe large. This was calculated and fitted using the numerical results after solving Eq. (61). We should note that when we are close to the critical Bingham number, yield stresses dominate over viscous stresses in the neighbourhood ofû = 0, hence the effect of μ( ) over B i,c is negligible. We have calculated B i,c for α μ = α B and α μ = 0 and we get the same surface as in Fig. 5a. Hence, for the rest of our work, we will concentrate on the effects on the variation of α B and fix α μ = 0.1.
The previous results can be used as a guide for the expected critical Bingham number for a given material. For example, using the physical parameters for lava flow reported in [9] we get α B /Pe ∼ O(1).

General topography, w s = 0
Now we are in a position to investigate the effects that topographical changes have on the spread of the material. The problem of gravity currents flowing on non-flat surfaces is a common occurrence and the ability to model and better understand the combined effects of viscoplasticity, topography, and heat transfer is important. These changes in geometry not only affect gravity currents or geophysical flows, but also for coating processes where the substrates are not flat but exhibit a topography that affect the thickness of the coating. Borrowing ideas from [27,41] we choose step-up and step-down topographies.
We use the following values for the constants in f (x) defined in Eq. (35): a = 0.1, b = 1.2 and c = 0.01. In [27], the authors showed that if the amplitude of the step is less than one, the difference in results for different magnitudes of the gradient of the topography (in our case O(1/c)) are negligible. By choosing a step with a small amplitude a we guarantee that the thin-film approximation for h remains valid. We will keep the isothermal Bingham number fixed at B i = 0.6. This value is close to B i,c for the temperature-independent rheology and we will consider changes in the ratio α B /Pe for the non-isothermal case.
In Fig. 6, we show profiles for h(t, x) and (t, x) for t = 1000 in a horizontal step-up topography. We consider α B /Pe = 0, 0.001, 0.01, 0.1, 1. As we can see, the only material that actually "climbs" the step is the isothermal one (α B /Pe = 0), see Fig. 6a. Thermal effects are considerable and even for the case α B /Pe = 0.001, the length of the block is noticeably affected by temperature, see Fig. 6c. In Fig. 6b we show the temperature profile for α B /Pe = 0.001, for all the other cases || || ∞ < 10 −10 . This was to be expected as we have shown that the ratio α B /Pe is the principal driver of the decay rate. Because of this rapid decay in temperature the material plateaus to its arrested state faster as α B /Pe is increased. Note also that for the non-isothermal cases, a yield point x y,∞ exists.
In Fig. 7 we repeat our simulations for the case of an inclined step-up topography. Results are similar to the horizontal case. Again, the isothermal case is the only one that climbs up the step. There are no noticeable  29) it is easy to see that the slope required for the whole block to deform should satisfy in steady state. For small values of m the effects on the length of block are marginal. In Figs. 8 and 9 we flip the topographical feature and investigate a step-down domain. Results are similar to the previous case with the main difference that the step down enhances the length of the isothermal block. For the horizontal topography the only non-isothermal case which climbs down the step is when α B /Pe = 0.001. As soon as thermal effects are strong enough, we have similar spreading rates as in the flat and step-up cases. When α B /Pe ≥ 0.01 a yield point x y,∞ exists.
We summarise our results in Table 1 and present steady state values of |x f,l | and x f,r − 1 for each value of the ratio α B /Pe. It is evident now why the steady state value of the length of the block remains the same for both horizontal and inclined step-up cases. As the right front moves forward due to the inclination of the topography, it "drags" the left front with it the same amount, the positive gradient of the topography and the yield stress inhibit the spreading of the material. This is obviously not the case for the step-down topography where gravity and topographical gradients enhance the movement of the right front of the block. As it has previously been considered in [23,29], we define a time-dependent fluid source, at constant unit temperature of the following form where H(x) is the Heaviside step function.
The main comparison will be done against the horizontal, flat case presented in Sect. 5.1. Even though we no longer have a steady state, we run our simulations up to t = 1000. In Fig. 10 we present the numerical results for height profiles and a time series of the block length. As in previous sections, we consider results for increasing ratio α B /Pe and fix B i = 0.6. The first thing to note is that the length of the block for all non-isothermal cases have plateaued to a fixed value. The material stops spreading laterally but not vertically. This is not the case for the isothermal case, because the Bingham number is below its critical value, viscous and gravity effects dominate the flow. As it is with the case for w s = 0, the material fully slumps and then spreads laterally helped by the source term, this is clearly seen in Fig. 10b in the black curve.
In Figs. 11a and b, we present the temperature profiles corresponding to the results shown in Fig. 10. In comparison with the results of previous sections the temperature behaves differently, even for α B /Pe = 1 the value of || || ∞ is non-negligible. Clearly, the inclusion of the source term (1 − )w s in Eq. (28) changes the rate of decay of the temperature from exponential to a power law, we confirm this in Fig. 11c where we plot || || ∞ as a function of time for α B /Pe = 1. As we can see In order to explore further the effects that the temperature-dependent rheology and source term have on the deformation of the material, we present snapshots of h(t, x) and h c (t, x) at times t = 0, 1, 10, 100, 1000 for α B /Pe = 0 in Figs. 12a and c, and for α B /Pe = 1 in Figs. 12b and d. For both cases and at early times, the material at the centre of the block yields and gets extruded upwards. In the isothermal case, at later times gravity effects dominate, the material is fully yielded and spreads laterally. No arrested state is reached and we have h c > 0 over the whole domain, see Fig. 12c. This behaviour was previously reported by Balmforth et. al. in [6] for isothermal lava domes.
The dynamics for the non-isothermal case are quite different. As we have shown before, the effective Bingham number for a temperature-dependent rheology is larger than B i . In this case, plastic effects dominate part of the flow. At early times the material gets extruded at the centre of the block and spreads laterally under gravity effects, see red curves in Figs. 12b and d. Quickly after, most of the block stops deforming, apart from the centre where the source keeps extruding material and in a small neighbourhood close to the fronts, (black curves). For later times, in the order of 10 2 time units and onward, the length of the block has reached a plateau and deformation is confined to the centre of the block.
In Fig. 13a, we show the time series h(t, 1/2). Clearly the peak grows but rather slowly. It seems that as t → ∞, because the material keeps being extruded, the size of the yield surface gets asymptotically smaller,   Fig. 13b if we compare the black curve (t = 1000) and the red curve (t = 10000). Intuitively, we can infer that the block will be fully yielded but it will take an infinite time to reach this state. It is worth noting that this behaviour is related to the temporal dependence of the source term w s (t, x). The choice of a power law behaviour is arbitrary. We could have chosen an extrusion rate which decays exponentially fast or terminates abruptly after a given time. This, of course, will modify the spreading dynamics of the block. These effects are being explored currently and will be presented in a future publication.

Conclusions
We have studied the flow of a block of a viscoplastic material with temperature-dependent rheology, spreading over different topographical configurations. The model used is based on the work of Bernabeu et al. [9] which consists of an evolution equation for the height of the material and a transport equation for the temperature. The model is valid for Pe > ε −1 and includes a Newton-type of cooling law that goes to zero asPe → ∞. We have assumed an exponential law for both the plastic viscosity and the yield stress as in [9]. We concentrate in a temperature-dependent yield stress and its effects on the spreading dynamics of the block of material. In Sect. 5.1.1 we studied how the critical Bingham number is affected by thermal effects. We developed an algebraic equation for an approximation of B i,c and through it we showed that the ratio α B /Pe would greatly affect the dynamics of the flow. Even for small values of the aforementioned ratio the spreading of the material is noticeably reduced. As expected, topographical features also perturb the spreading of the material but unless the slope of the inclined plane or the size of the step are large enough to counter/aid the effects of gravity, yield stress and thermal effects dominate the flow. In Sect. 5.3 we consider the case when we have a time dependent source of fluid for which its intensity decreases as time increases, i.e. w s → 0 as t → ∞. In this case, thermal effects are of utter importance, if we do not consider them and the Bingham number is below its critical value the block does not arrest. On the other hand, if the yield stress is temperature-dependent, its effective Bingham number will be above B i,c and the deformation of the material is practically confined to the centre of the domain. To our knowledge, this is one of the first attempts to characterise the flow of a block of viscoplastic material driven by gravity where topographical and non-isothermal rheology are considered. The next step in our investigations will be to compare numerical simulations with experiments. Experimental data for nonisothermal gravity-driven flows is available, to our knowledge, only for Newtonian fluids. For viscoplastic materials, only flows under isothermal conditions have been studied, see [13] for an extensive review of benchmark problems for lava-flow models. Bernabeu et al. [9] validated their model and numerical results against the experimental data (Newtonian) of Garel et al. reported in [18] and agreement is good. Nonisothermal experiments using a viscoplastic material are being performed by members of our group at the University of Canterbury. Our final goal is twofold: to study how thermal effects affect the spreading of yieldstress materials; and given surface data, predict rheological properties of the fluid. As we have mentioned earlier in this paper, lava is known to have a temperature-dependent rheology and can be the focus of a particular application in a future work. Rheological parameters reported in [9], show that the ratio α B /Pe ∼ O(1) for some lava flows. As shown in the present work, for small values of the aforementioned ratio, non-isothermal effects are non-negligible. Many questions remain and we hope to answer them in future publications.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions This work is part of the project "Indirect measurement of lava rheology", Marsden Fund UOC1802. The funding is kindly acknowledged.

Declarations
Conflict of interest This work has no competing interests that might be perceived to influence the results and/or discussion.

Ethical approval Not applicable.
Data availability For any access to data or material presented in the manuscript, contact the corresponding author.
Author contribution MM-G wrote the main manuscript text and all the authors reviewed and discussed the manuscript.