Far-field boundary conditions in channeled lava flow with viscous dissipation

10 Cooling and dynamics of lava flowing in a rectangular channel driven by the gravity force is 11 numerically modeled. The purpose is to evaluate the thermal process as a function of time involving 12 the liquid lava in contact with the solid boundary that flanks lava. Lava rheology is dependent on 13 temperature and strain rate according to a power law function. The model couples dynamics and 14 thermodynamics inside the lava channel and describes the thermal evolution of the solid boundary 15 enclosing the channel. Numerical tests indicate that the solution of the thermo-dynamical problem 16 is independent of the mesh. The boundary condition at the ground and at the levees is treated 17 assuming a solid boundary around the lava flow across which lava can exchange heat by 18 conduction. A far field thermal boundary condition allows to overcome the assumption of constant 19 temperature or constant heat flow as boundary conditions, providing more realistic results. The 20 effect of viscous heating is evaluated and discussed. 21 22

When the problem of lava flowing under the effect of gravity is resolved, the power law rheology introduces a non-linearity in the diffusion term of the momentum equation and consequently an analytical solution of the differential equations governing the motion does not exist.Furthermore, if the viscosity function also includes temperature dependence, the thermal equation is coupled to the dynamic equation.Hence the need to solve nu-merically the thermo-dynamic equations describing the flows of a fluid as complex as lava.
In the thermal modeling of lava flows, the greatest difficulties are due to the different thermal exchanges both external (surface thermal radiation, forced convection, conduction at the base) and internal (axial advection, viscous dissipation, latent heat, internal conduction) to be taken into account.Despite the great progress made in numerical modeling, it is necessary to assume some simplifying hypotheses.
The need for simplifying assumptions in numerical models of lava flows is described in review articles [Costa and Macedonio, 2005b;Cordonnier et al., 2015;Dietterich et al., 2017].The authors examine the most relevant works on the numerical modeling of lava flows and what emerges is that, due to the high complexity of transport equations, the numerical solution of the complete threedimensional problem for real lava flows is often in-tractable.This has led to concentrate the main efforts on the development of software able to quickly describe the evolution of a lava flow, based on simplified theoretical models, for the purpose of volcanic hazard assessment.
Parallel to the development of general numerical models closely linked to volcanic hazard, it is necessary to study the details of what happens when, during an effusive process, a fluid like lava flows in a channel.Filippucci et al. [2017] discussed the problem of viscous dissipation in channeled lava flows and found that the effect of viscous heating is irrelevant for the most of the channel domain apart from the boundaries where some temperature differences where noticed.These temperature differences brought the authors to hypothesize that the fluid motion can lose the laminarity at the boundary (since locally the Reynolds number exceeds the laminar/turbulent threshold).The same effect was found by Costa and Macedonio [2005] who found that the fluid can develop secondary flows at the boundaries as an effect of viscous heating.
In dealing with problems involving the heat equation, in this case coupled to the momentum equation by the rheological function, the boundary condition for a channeled flow is classically chosen between constant temperature and constant heat flow [Patankar, 1980;Ferziger and Peric, 2002].
From the observation that the effect of viscous heating is observed only at the boundaries, in this paper we considered a third possibility, that is a conductive solid boundary with which the fluid can exchange heat.In particular, in this paper we focus on the thermal process involving the liquid lava in contact with the solid edges of the channel, assuming that the fluid can cool down by losing heat by conduction in the hosting rocks and by radiation in the atmosphere and can heat up by the effect of heat advection and viscous dissipation.Some authors [Costa et al., 2007;Filippucci, 2018] used far-field boundary conditions to model the thermal interaction between the hot fluid flow and the host rock with numerical methods.Similarly, in this work, the solid boundary condition is treated realistically considering that the rock that hosts the lava flow can exchange heat by conduction.

MATHEMATICAL PROBLEM
The constitutive equation for a power law fluid is (1) where σ ij is the stress tensor, e ij is the strain rate tensor, k is the fluid consistency, n is the power law exponent which is a measure of nonlinearity, and (2) where I 2 is the second invariant of the strain rate tensor.The apparent viscosity η a of the fluid is (3) If n is lower than 1, the fluid is pseudoplastic and it thins with an increase in stress.If n is equal to 1, the fluid is Newtonian.If n is greater than 1, the fluid is dilatant and it thickens with an increase in stress [White, 2005].
We assume a viscous fluid flowing in the x direction in a rectangular conduit inclined at an angle α, perpendicularly to the section of the conduit in the yz plane.The channel is surrounded by edges of solid material of thickness a s with which it can exchange heat by conduction.
The channel width is 2a l , thickness h l and length L. The sketch of the model with the coordinate system is shown in Figure 1 and the values of the parameters are listed in Table 1.
The flow is laminar and subjected to the gravity force.Pressure changes are negligible with respect to body forces.The velocity is approximately axial and varies with the lateral coordinates v x (y,z), v y =0, v z =0.The fluid is isotropic and incompressible with constant density ρ, thermal conductivity K, specific heat capacity c p .
The equation of motion in the transient state for a gravity driven flow down an inclined rectangular channel is [Filippucci et al., 2013a]: where v x is the x component of velocity, g is the gravity acceleration and α is the slope angle.The apparent viscosity is [Filippucci et al., 2005]: (5) where v x is the x component of velocity and both the fluid consistency k and the power law exponent n depend on temperature.The temperature dependence of k and n is given by Hobiger et al. ( 2011): (6) (7) where k 0 , p 1 , p 2 , p 3 and p 4 are constant parameters listed in Table 1.
The numerical problem of the dynamics of lava flows was already detailed, by studying the sensitivity of the solution to the variation of the power law exponent n and to the fluid consistency k [Filippucci et al., 2010;2011;2013b], by using a temperature-dependent power-law rheology and analyzing the thermal effects due to heat advection [Filippucci et al., 2013a] and to viscous dissipation [Filippucci et al., 2017].The problem of the solid boundaries has been considered by Filippucci [2018].
We assume that the fluid flow is transient, laminar and subjected to the gravity force.Downflow pressure changes are negligible with respect to the body force.Since the channel cross section is rectangular and the Reynolds number is low almost everywhere in the domain [Filippucci et al., 2017], the assumption of laminar flow implies that the velocity is approximately axial, but it may depend on all the coordinates: v x = v x (x,y,z).Differently to the approach of Filippucci [2018], we include the effect of viscous dissipation as an internal heat source.We neglect the effect of the latent heat of crystallization/fusion.
The boundary conditions are the no-slip at the walls and the zero-stress at the top of the flow.Since the solution is computed in a half domain of width a/2, thickness h and length L, we consider as boundary condition the symmetry of the problem with respect to the xz plane.So, the boundary conditions are the following: T h e initial condition for velocity is the steady state numerical solution at the initial effusion temperature v x (x,y,z, t=0) = v x (T e ) .
The time dependent heat equation inside the channel takes into account the effect of thermal exchange by heat advection, conduction and viscous dissipation: We can neglect the effect of thermal conduction in the flow direction as it is of secondary importance with respect to thermal advection in the flow direction [Filippucci et al., 2013a].
The time dependent heat equation outside the channel, in the solid boundary, is purely conductive: The thermal boundary conditions are the assumptions of a radiative heat flux q r at the upper flow surface, a constant temperature T w at the outer solid walls (such a far field condition is imposed at a distance equal to three times the half-width of the channel), a constant effusion temperature T e at the vent and the symmetry of the problem with respect to the xz plane: where q r = σεT u 4 and σ is the Stefan-Boltzmann constant, ε is the surface emissivity of lava and T u is the temperature of the upper surface at z = 0 (the atmospheric temperature is assumed negligible with respect to T u ).
At time t = 0, the liquid lava has a uniform temperature T e , the velocity is the stationary solution of the dynamic equation with T = T e and the outer solid boundary has uniform temperature T w .The choice of the initial condition was made in order to compare this solution with that of Filippucci et al. (2017).Moreover, the assumption of the extrusion temperature of the fluid as starting condition in numerical studies is widely used and accepted in finite element/volume modeling [Costa and Macedonio, 2003;Patrick, 2004;Bernabeu et al., 2016;among others].
At time t > 0, the radiative heat flux q r , the far field constant temperature T w and the constant effusive temperature T e are imposed.Since x = L is the outflow boundary and both temperature and velocity need to be computed there, no boundary condition at x = L is necessary.
The dynamic and the thermal equation are coupled by the temperature dependence of viscosity in the dynamic equation and by the viscous dissipation term in the heat equation.The algorithm is written by the authors in Fortran language.The space discretization is obtained by the control volume integration method [Patankar, 1980] using a static mesh approach and power law interpolation functions between the nearest grid points.The radiative condition at the boundary z = 0 depends on the fourth power of temperature and needs to be linearized in order to be treated as a source term of the heat equation.For the discretization of the transient term, the integration over the time interval t is made by using a fully implicit scheme.The iterative solutor of the discretized equation is the Gauss-Seidel one [Filippucci et al., 2013a, for details of the flow chart procedure].
The solution is tested to verify the independence of the mesh and the test is shown in Figure 2. The computational problem was solved by considering three grids of different sizes (52×52, 102×102, 202×202) to discretize the (y, z) section, transversal to the fluid direction.The space discretization along the x direction was fixed to 101 control volumes.The time solution is stopped at t = 10 6 s since for long times the temperature approaches the steady-state solution.The channel geometry for the numerical test is described by the parameters in Table 1.The temperature profile along the z coordinate at the outflow boundary slightly varies with varying the mesh size, indicating that the problem is independent of the control volume size.
As expected, the finest mesh (201×201) needs a very large computational cost to achieve the convergence, which means a long time for calculation.In the following, for the problem with geometrical parameters listed in Table 1, the mesh y × z × x = 102 × 102 × 101 was used as the best compromise between accuracy and computation time.In particular, in this case, the computational time is approximately 2 days of computation for the problem with the viscous dissipation term, and approximately 1 day for the solution without the viscous dissipation term.

RESULTS AND DISCUSSION
The problem is illustrated in Figure 1 with the physical and geometrical parameters given in Table 1.We have evaluated the temperature distribution in an inclined lava channel flanked by solid levees, with which the flowing lava interacts thermally: while lava cools down, the levees are heated in turn.Results are obtained considering that the lava cools by thermal radiation at the free surface and by thermal conduction at the solid boundaries.Moreover, the lava can be heated by the effect of viscous dissipation.The solid surrounding rocks, in turn, heat up by thermal conduction.
The problem is transient, and the first 10 6 s (approximately 1 day) are modeled, assuming that at t = 0 the whole channel has a temperature T equal to the effusion temperature T e and the levees have uniform wall temperature T w .
The choice of the distance of the constant temperature boundary condition is arbitrary.Costa et al. [2007] used a far-field temperature at an arbitrary distance of 10 times the radius of the magma conduit.Our choice is dictated as a compromise between far-field and computational time costs.This choice can be considered adequate since, after t = 10 6 s of cooling, the boundary is still not heated, as it can be seen in Figure 3a.
In a previous paper, Filippucci [2018] has numerically solved the same problem, but neglecting the viscous dissipation term in the heat equation.In Figure 3b, we plotted the difference with the solution of Filippucci (2018) in terms of temperature T (y, z) map in the vertical cross section at the outflow boundary (x = L) at t = 10 6 s.It can be observed that the viscous dissipation term has an effect in the parts of the lava channel in contact with the boundaries and the ground.If we consider the whole lava channel, the effect of viscous heating can be considered of secondary importance with respect to the advective term.If we consider the behavior of the fluid in contact with the solid edges, the heat addition due to viscous dissipation can bring the fluid flow to change the motion from laminar to locally turbulent [Filippucci et al., 2017] and to develop secondary flows [Costa andMacedonio, 2003, 2005a] In Figure 4, we plotted the temperature variation with time at 8 monitoring points (P1, ..., P8) selected inside the lava channel and corresponding to the black points in Figure 4A.If we observe the temperature evolution with time at some points of the channel section in the flow outlet area (Figure 4), we realize that there are fluid regions that are almost at constant temperatures equal to those of T e and areas with temperature that oscillates without differences between dissipative and non-dissipative case with the exception of point P5 in contact with the ground at the center of the channel.At points P1, P6 and P8, temperature remains constant and equal to that of effusion T(t) = T e during the prescribed time.At these points, the temperature difference between the dissipative and non-dissipative case are in the order of tenths of a centigrade degree.At point P2, temperature decreases monotonically during the first 10 3 s approximately and then remains constant for the following time.The effect of viscous dissipation is to dampen the cooling process and to allow the fluid to be at a constant higher temperature in less time.At points P3, P4 and P7, the temperature has an oscillating behavior in time: at first it decreases and then returns to increase without significant differences between the dissipative and nondissipative case and with a greater oscillation amplitude at point P4 than in the other points.At point P5, the temperature behavior with time is similar to the one just described, but in this case the difference between dissipative and non-dissipative case is sharper.
The main difference with the work of Filippucci et al. [2017] and then between boundary conditions with constant temperature gradient and far field boundary conditions, where temperature gradient is free to vary as a function of the thermal conduction rate, is just shown in Figure 4 (points P3, P4, P5 and P7).
In fact, at the same points, in the case of boundary condition with constant heat flux, temperature increases monotonically in the presence of viscous dissipation, while it decreases monotonically in the absence of dissipation [Figure 7 in Filippucci et al., 2017].Considering far field boundary conditions, the temperature initially decreases and after a certain time internal it starts to increase and this oscillation is independent from the viscous dissipation term in the heat equation.
This oscillation can be interpreted considering that, as the lava flow is emplaced, the difference in temperature between the host rocks and the hosted lava is so high that it causes a very intense conductive heat flow.This flow cools the lava in contact with rocks, as shown by the descending part of the curves in Figure 4 (points P3, P4, P5 and P7).Over time, the host rocks heat up as a result of heat conduction and the heat flow begins to decrease and so does the cooling rate of the lava in contact with rocks, leading to the minimum of the curves in Figure 4 (points P3, P4, P5 and P7).Continuing with time, the conductive heat flow becomes no longer the dominant mechanism, since advective heating begins to prevail due to the isothermal core at the center of the channel, that flows with temperature equal to that of effusion.This change in heat transfer mechanism leads the hosted lava in contact with the host rocks to heat up, as shown by the ascending part of the curves in Figure 4 (points P3, P4, P5 and P7).The effect of viscous dissipation during the described thermal process is that of reducing the oscillation amplitude Figure 4 (points P5).

CONCLUSIONS
The purpose of this work is to investigate how the choice of boundary conditions may affect the results of a model of a channeled lava flow, in particular by evaluating the importance of the viscous dissipation term.To do this, we developed a flow model with far-field constant temperature boundary conditions, allowing the liquid lava to exchange heat with the host solid rock.The results were compared with those obtained in a previous work [Filippucci et al., 2018] in which the cooling at the levees and at the ground was modeled by imposing a constant conductive heat flow.With the model presented in this paper, we go beyond the classical boundary conditions at the channel walls, which assume constant temperature or constant heat flow, usually adopted in works dealing with lava flow simulation, as reviewed by Costa and Macedonio [2005b].Far field boundary conditions are more realistic since it is not necessary to impose a constant temperature nor an arbitrary constant temperature gradient at the channel levees and at the ground.
We solved numerically the dynamic and heat equations of a lava flowing by gravity inside a channel with rectangular cross section, flanked by a thick solid levee, by using the finite volume method [Patankar, 1980].
The dynamic equation is solved only in the liquid domain, while the heat equation is solved both in the liquid and in the solid domain.The viscous dissipation term is considered in the heat equation.The solution was tested in order to verify that the convergence of the numerical problem is independent of the mesh size.The results indicate that the solid edge interacts with the liquid lava, both at the levees and the ground, causing an initial cooling due to heat conduction and a subsequent heating due to heat advection.This thermal process is not affected by viscous dissipation, which only acts by decreasing the temperature variation interval between the cooling and the heating process.On the contrary, when the boundary condition at the levees and at the ground is a constant temperature gradient [Filippucci et al., 2017] or a constant temperature [Costa et al., 2007], the effect of viscous dissipation at the boundaries appears to be of great importance, causing an increase of the Reynolds number from laminar to turbulent values [Filippucci et al., 2017] and triggering secondary flows [Costa et al., 2007].
Filippucci [2018], using the same far field boundary conditions but neglecting viscous heating, found that the channel levees can melt, since they can heat over the solidus temperature.In this work, including viscous dissipation, we observe the same thermal behavior of the edge (Figure 3a), but the effect of thermal erosion cannot be evaluated quantitatively.In order to analyze melting of the solid edges and changes in the channel morphology, we should adopt a moving boundary, so as to account that portions of the host rocks may pass from solid to liquid behavior and, vice versa, portions of lava in the channel may pass from liquid to solid behavior.As reviewed by Costa and Macedonio [2005b], works dealing with lava flows impose that the channel boundaries are taken at constant temperature or at constant temperature gradient.Differently, far field thermal boundary condition allows to overcome this simplification and makes the physical treatment of the problem more realistic.

3FARFIGURE 1 .
FIGURE 1. Coordinate system and geometrical parameters.Boundary conditions are also indicated.

FIGURE 2 .
FIGURE 2. Temperature profile along z coordinate at x=L, y=0, t=10 5 s, for different mesh sizes, as indicated in the figure legend.

FIGURE 4 .
FIGURE 4. Temperature vs time at different monitoring points (P1,..,P8) of the lava channel, as indicated in (A).The label on each plot corresponds to the point in (A).Solid line: solution of the heat equation considering the viscous dissipation term; dashed line: solution of the heat equation neglecting the viscous dissipation term [from Filippucci, 2018].

TABLE 1 .
Values of the fixed model parameters