Numerical simulation of lava flow using a GPU SPH model

A smoothed particle hydrodynamics (SPH) method for lava-flow modeling was implemented on a graphical processing unit (GPU) using the compute unified device architecture (CUDA) developed by NVIDIA. This resulted in speed-ups of up to two orders of magnitude. The three-dimensional model can simulate lava flow on a real topography with free-surface, nonNewtonian fluids, and with phase change. The entire SPH code has three main components, neighbor list construction, force computation, and integration of the equation of motion, and it is computed on the GPU, fully exploiting the computational power. The simulation speed achieved is one to two orders of magnitude faster than the equivalent central processing unit (CPU) code. This GPU implementation of SPH allows high resolution SPH modeling in hours and days, rather than in weeks and months, on inexpensive and readily available hardware.


Introduction
The development of physical-mathematical models that can describe the spatial and temporal evolution of natural phenomena is a key point in many scientific areas.The ability to predict potentially affected areas of high-risk from volcanic phenomena, such as lava flow, is essential to support risk mitigation and land planning, in combination with laboratory and field observations.Several physical models and numerical methods have already been applied to simulate lava flow under some simplified assumptions.Models based on the concept of maximum slope and stochastic perturbation of topography [Favalli et al. 2005] have been integrated with models that include a more complete physical description.These have been implemented with cellular nonlinear networks [Del Negro et al. 2005] and cellular automata, which can describe the spatial and temporal evolution of lava flow on the basis of given eruptive parameters [Vicari et al. 2007, Del Negro et al. 2008, Hérault et al. 2009].
Over recent years, these models have been applied with success in collaboration with the Department of Civil Protection for the creation of hazard scenarios during Mount Etna eruptions.However, they are inadequate for the description of more sophisticated phenomena, like crust and lava-tube formation, and ephemeral vent opening.The occurrence of these phenomena can strongly increase the hazard associated with lava flow, as seen, for example, during the 1669 Etna eruption in which the city of Catania (Italy) was destroyed, and in the 1981 and 1991-1993 Etna eruptions when lava threatened the towns of Randazzo (Italy) and Zafferana Etnea (Italy), respectively.For the purposes of the Department of Civil Protection, this aspect has significantly grown in importance with the urban expansion that has almost entirely covered the Etnean foothill area more recently, which now provides permanent residence to about 800,000 people.Even though they do not represent an immediate threat to people's lives, Etnean eruptions are a constant menace to housing and infrastructure.The scientific community must address this problem by the provision of working models of the volcano that have improved precision, so these can help determine the areas of Mount Etna that are exposed to the greatest risks of lava invasion.Thus, knowledge of the processes for collocation and expansion of lava flow is a fundamental requisite for the evaluation of lava-invasion risk.
A complete simulation of lava flow is challenging from a modeling, numerical and computational point of view.The natural topography irregularities, the dynamic free boundaries, and phenomena such as solidification and friction, presence of floating solid bodies or other obstacles, and their eventual fragmentation, all make this problem difficult to solve using traditional numerical methods, such us finite volumes or finite elements.Application of these approaches to simplified cases can be found in Quareni et al. [2004], Hale [2008] and Filippucci et al. [2010].Eulerian discretization has an additional problem with boundary tracking, which for complex fluids such as lava, is further complicated by the presence of the internal boundaries from the solidification fronts.Lagrangian methods, such as finite elements, are instead challenged by problems that relate to the highly deformable nature of lava flow, which inevitably leads to significant deformation in the finite element structure, with the consequent loss of accuracy and the need for periodic remeshing.
A new class of methods, which are known as meshless, mesh-free or particle methods, has emerged over the past few decades as an alternative to traditional grid-based methods.Smoothed particles hydrodynamics (SPH) was developed in the 1970s by Lucy and by Gingold and Monaghan for astrophysical simulations [Lucy 1977, Gingold andMonaghan 1977].Then, it was extended to fluid dynamics problems [Monaghan 1992], such as free-surface problems [Monaghan 1994] and viscous flow [Morris et al. 1997].Moreover, SPH was applied to thermal conduction problems [Cleary and Monaghan 1999] and phase transition [Monaghan et al. 2005].These belong to the class of Lagrangian meshless methods in which the problem to be solved is discretized using particles that are free to move, rather than using fixed grids or meshes.The governing partial differential equations are converted into ordinary differential equations for the evolution of the particle parameters (velocity, position, temperature, among others).Many of the advantages of particle methods are of direct application to lava-flow modeling; e.g. the ability to simulate free surface flow, multi-phase physics, and solid/fluid interactions.Another unique and attractive characteristic of particle methods is that particles can carry material properties, rather than acting only as interpolation points, which is instead the case for other meshless methods.
The major obstacle for using particle methods in modeling natural phenomena is the high computational resources required, which is traditionally solved by using expensive clusters.On the other hand, these methods typically show a very high degree of parallelism, which allows their implementation on parallel computing hardware to be very efficient.To achieve high computing performance, we used the low-cost, energy-effective parallel-computing capabilities offered by the new generations of graphic processing units (GPUs).A single GPU contains a large number of computing cores that can access a global memory, and it is designed to process large quantities of data in parallel.These characteristics make GPUs particularly appropriate for the implementation of highly parallelizable algorithms (such as those needed for particle methods), and allow them to achieve performances of up to two orders of magnitude higher than those of standard CPUs.
In the present study, we show a recent application of SPH to the modeling of lava flow, which was developed at the Istituto Nazionale di Geofisica e Vulcanologia (INGV), Sezione di Catania, within the course of the LAVA Project.
This three-dimensional model can describe the flow of a fluid with thermal-dependent nonNewtonian behavior, which includes phase transition, on a natural topography and with a free surface.This allows for the accurate forecasting of the possible paths of a lava flow during an eruption.

Notation
Throughout this manuscript we use the following conventions: • Vectors will be written in bold; • The reference coordinate system is cartesian orthogonal with axes x 1 , x 2 and x 3 ;

Mathematical principles of smoothed particles hydrodynamics
Within the SPH method, the fluid is represented as a collection of "particles", whereby each particle represents a portion of the fluid.The motion of the particles is driven by forces that are computed following the Navier-Stokes equations, which can be coupled with other equations to describe, e.g., the evolution of temperature.
First, let us consider a scalar field f defined in a domain X.We can write: (1) Let us give an approximation of d with function W with compact support 1 , where the diameter is kh, satisfying: (2) (3) The function W is known as the kernel, and h is the smoothing length.By substituting Equation (3) in Equation (1), we obtain: The previous integral, denoted as Gf (r)H, is the SPH approximation of f(r) with the kernel W: (4) Secondly, we can compute an approximate value of Gf(r)H using a discrete set of particles that are positioned at r 1 , …, r N and that occupy the volumes V a , respectively, and for which the value of the function f is known, as: Basically, the sum includes only the particles lying inside the support of W. Moreover the volume V a occupied by the particle a is equal to the ratio between its mass, m a , and density, t a .We obtain: (5) From this point on will be substituted by =, and For the density, there is, for instance: (6) In brief, the SPH method is build on two successive approximations: • The kernel approximation, in which we replace the Dirac distribution in Equation (1) by the kernel W leading to Equation (4); • The particle approximation, in which we approximate the integral of Equation ( 4) by a discrete sum over the particles, leading to expression of Equation ( 5).
The expressions of Equations ( 1) and ( 4) are convolution products between f and a regularization function, and therefore they have the property of convolution products.In particular, the convolution product of the derivative of f with the regularization function is equal to the convolution product of f with the derivative of the regularization function: Thus, with the SPH notation, we obtain: (7) This relation is valid if the integral on ^X is zero, and so only if the support of W is included in X.This condition is not verified for particles near solid boundaries or near the free surface.This is one of the main difficulties encountered for boundary treatment with SPH.To overcome this, for the particles near solid boundaries we can: • add a corrective term that represents the integral on ^X, which in this case is not zero; • add some particles outside of the solid boundaries, to fulfill the conditions on the support of W.
The function W must satisfy the conditions of Equations ( 2) and (3), and moreover, Equation (6) implies that W is positive.In practice, the kernels are also chosen symmetrically: W(r − r´,h) = W(||r − r´||,h).A detailed description of how to build kernel functions can be found in Liu and Liu [2003].
It is possible to obtain different SPH formulations, with the ones that operate symmetrically on the particles usually preferred, because these tend to improve the preservation of the laws of conservation.For example, in the case of the gradient, we can write: (8) Applying the SPH method to a partial differential equation, for every particle we obtain an ordinary differential equation.For example, using Equation (8), we have the widely used SPH form of the continuity shown in Equation ( 9): The main advantages of SPH are: • It is a purely Lagrangian method.• There are no limitations prescribed by problem geometry.
• There are no limitations prescribed by the differences between the final and initial states.
• The particle density is linked to actual fluid density, and this means that the precision increases in the area with the larger density.
The main drawbacks of SPH are: • A special treatment of solid boundaries is required: as we cannot directly solve the incompressibility condition within the SPH framework, it is necessary to introduce a state equation and to treat an incompressible fluid as a quasicompressible one, or to solve independently.• There is decreasing precision in the area with the lower density.

The equation of the problem
In this report, we choose to use the SPH without independently solving .So, even if lava is incompressible, this choice will lead to treatment of the lava as a quasi-compressible fluid, by the introduction of a state equation (see Section 4.3).The liquid phase of the lava is then modeled by the Navier-Stokes equations (continuity equation, forces balance): (9) (10) coupled with the heat equation: With a minimum viscosity of 0.1 m 2 s −1 , a maximum velocity in the order of 1 ms −1 and a hydraulic radius of the order of 50 m, we obtain a Reynolds number of 500.Of course, the previous assumptions on viscosity, velocity and hydraulic radius are a worst-case scenario, and the given Reynolds number is an absolute maximum.For example, in the simulations that we will show for the numerical experiments, the minimum viscosity is 1.2 m 2 s −1 , so with the same assumptions on velocity and hydraulic radius, we get a Reynolds number of less than 50.We can therefore assume that the Etna lava flows are laminar.

Rheological model -constitutive equation
The constitutive equation is the relation between the viscous stress tensor x ij and the strain tensor e ij .In the case of lava flow, the parameters of the constitutive equation depend on the temperature.
The choice of the rheological model and the properties of lava flow are a nontrivial question.To address this topic, different rheologies have been taken into account in our SPH lava-flow model: a Newtonian model and a Bingham model [Pinkerton and Norton 1995, Giordano and Dingwell 2003, James et al. 2004], where both are dependent on the temperature and the water content of the lava.Moreover, we also take this opportunity to introduce a general model that encompasses many nonNewtonian rheologies, known as the Herschel-Bulkley model [Balmforth et al. 2000].

The Herschel-Bulkley model
The constitutive equations for a Herschel-Bulkley fluid are: (12) where k is the consistency index, n is the power law index, x 0 is the yield stress, and A II is the second invariant 2 of the tensor A. In practice, we have a threshold fluid with a solid behavior under certain stress conditions.
In our SPH model we use the regularized model proposed by Zhu et al. [2005], which is an extension of the Papanastasiou [1987] model for Bingham fluids: with: (13) where parameters n 0 , t 1 and m are measured in Pa • s, s and s, respectively, and are chosen to obtain a large apparent viscosity for a vanishing shear stress.The main advantage of this model is to eliminate the discontinuities shown in Equation ( 12) and the associated numerical difficulties.From this point, this model is referred to as the modified Herschel-Bulkley model.
] ] The Bingham model A Bingham fluid is a Herschel-Bulkley fluid with a power law index n = 1.So, the constitutive equations are Equation ( 12) with n = 1.The corresponding regularized model from Zhu et al. [2005] is Equation (13) with t 1 = 0. Obviously, we recover the original Papansatasiou model: (14) From this point on, this model is referred to as the modified Bingham model.

The momentum equation used in the model
Treating nonNewtonian rheologies in the previous way, the momentum Equation (10) becomes: (15) This equation is formally equivalent to that obtained for a Newtonian fluid.The only difference is that the apparent viscosity n app is computed as shown previously.

The thermal aspects
We must take into account the phase transition and the two boundary conditions of: • Radiative losses at the free surface, which are given by: (16) • Conduction to the ground.As a first approximation, convective losses at the free surface can be neglected with respect to radiative losses [Harris and Rowland 2001], and we can neglect heat diffusion into the ground i.e the ground temperature remains constant.Phase transition will be treated without the introduction of any evolution equation of solidification fronts, as explained in Section 4.3.

SPH formulation of the Navier-Stokes equation
For the continuity equation, using the gradient expression (8) leads to: For the pressure part of the momentum Equation (15), we use the SPH formulation reported in Monaghan [1992], and for the viscous part, we use the equation derived by Morris et al. [1997]: Obviously, n a is the apparent viscosity of particle a computed in agreement with the chosen rheological model.For a nonNewtonian model, the apparent viscosity depends on the second invariant of the stress tensor, which was previously computed using Equation ( 8), to retrieve an SPH approximation of the velocity gradient.
As we are considering the fluid as quasi-compressible, we need to introduce a state equation to close our model.We use the Tait equation [Monaghan 1994], which is commonly used for free surface problems: with c = 7. Parameter B is chosen to guarantee that the speed of sound c 0 for the density at rest t 0 is 10 times greater than the maximum velocity of the flow.Here, we obtain:

SPH formulation of thermal aspects
For the conduction part of the heat Equation ( 11), we use the SPH formulation proposed by Cleary and Monaghan [1999], and for the viscous heating part, we use the SPH formulation found in Cleary and Ha [2003]: (19) where p = 4.963.
To treat phase transition, we use a slightly modified version of the method proposed by Monaghan et al. [2005]: • We add to each particle the parameter q, which represents the fraction of latent heat released during phase transition (if q = 1, the particle is solid; if q = 0, the particle is liquid).
• During the evolution of the temperature of a particle computed with Equation ( 19), if we have T > T s with q < 1, q is increased by c p (T -T s )/L, and T remains fixed to T s .
• If q ≥ 1, the particle becomes solid, and its temperature is given the value , after which it evolves normally according to Equation (19).
In this way, the different solidification fronts are formed implicitly, and they evolve automatically.These are represented by the particles with 0 < q < 1.Here, T s is the solidification temperature.However we cannot define the phase change temperature for lava, so instead we use the solidus temperature.

Ground interaction and boundary conditions
A way to treat solid boundaries is to fill them with particles that only generate a repulsive force [Monaghan and Kos 1999]: where r 0 is the radius of influence of the repulsive force (typically equal to the initial inter-particle distance Dp), p 1 = 12, p 2 = 6 and D is a problem-dependent parameter.
In the case of flow on a real topography, this approach is not convenient: we would have to generate a large number of particles to cover the area that would potentially be invaded by lava flow, even though only a very small number of these will actually interact with fluid particles with force, as Equation ( 20).Therefore, we have implemented a new method to compute a ground normal repulsive force (Figure 1): • Given a point P(x, y, z), we use linear interpolation to compute the heights z 0 , z 1 and z 2 that correspond to the points (x, y), (x + Dx, y) and (x, y + Dy).
• The ground normal vector n is taken as equal to the normal vector to the plane defined by the three points (x, y, z 0 ), (x + Dx, y, z 1 ) and (x, y + Dy, z 2 ).
• For every particle near the ground, we apply a repulsive force ||f||n.
In this way, starting from a digital elevation model (DEM) that even has sparse resolution, we can generate regular normals with a higher resolution (Figure 2).
The classic boundary particle approach has two significant limitations.First, the repulsive force depends not only on the distance to the boundary, but also on the boundary particle distribution.Secondly, the repulsive force is not exactly normal to the boundary.With our approach, we obtain a strictly normal force that only depends on the distance of the particle to the boundary.

Ground friction
Lava is a highly viscous fluid, so we have to introduce the effects of ground friction.Starting from the definition of viscosity, we apply the following to every particle in contact with the ground a force: (21 where u T is the velocity component that is tangent to the ground (Figure 3), and S 1 is the contact surface.
A particle is considered to be in contact with the ground when its distance from the ground is less than r 0 .A comparative study between the simulations of a viscous flow on a slope and the corresponding analytical solution leads us to choose: where S g is the surface of the portion of the flow in contact with the ground, and N g is the number of particles in contact with the ground.The evaluation of S g is explained in Section 5.1.

. Free surface detection
To deal with radiative losses at the free surface, we need to detect the particles that belong to the free surface.In our case, the lava flow is a low Reynolds number flow and has a very regular free surface.Taking into account these features we can detect free surface particles with the following algorithms: • For every particle a we consider C a as the cone with vertical axis and vertex on particle a; let i be the angular aperture of the cone.
• If C a does not contain any particles, the particle a is on the free surface (Figure 4, particle a).
• Otherwise a does not belong to the free surface (Figure 4, particle b).

Boundary conditions for the heat equation
We have two boundary conditions for the heat equation: • Radiative losses at the free surface −ln.∇T = = k b f • Heat flux continuity across the lava/terrain interface.Theoretically these two conditions result in a nonlinear system and a linear system, respectively, that need to be solved at each time step along with the solution of the heat equation into the ground.Practically, as the time step fulfilling the global Courant-Friedrich-Lewy (CFL) condition is a few orders of magnitude lower than that required by the thermal diffusion part (as explained in Section 4.5), and as we consider that the contact temperature of lava to the ground remains constant, we simply consider that during Dt: • a free surface particle exchanges a heat quantity of: which results in a temperature change of: where S 2 is the surface participating in the radiative transfer.
In a similar way to what was done for the ground friction, we state that S 2 = S f /N f where N f is the number of particles belonging to the free surface and S f its extent.
• a ground particle exchanges a heat quantity equal to , which results in a temperature change equal to: (24 where r is the distance of the particle from the ground.

Numerical integration
We use an explicit second-order predictor corrector scheme.To ensure its numerical stability, the time step has to fulfill the CFL condition, which in our case is made up of four parts.Within a time step Dt: • a force on a particle must not induce a displacement greater than a fraction of the smoothing length ; • a signal must not propagate more than a fraction of the smoothing length ; • viscous diffusion must not propagate more than a fraction of the smoothing length ; • heat diffusion must not propagate more than a fraction of the smoothing length .So the global CFL condition for this problem is: For all of the simulations we use C 1 = C 2 = 0.3, C 3 = 0.125 and C 4 = 0.1.For the typical values of the parameters of Etnean lava (Table 1) and as the maximum flow velocity is of the order of 1 ms −1 , the thermal part of the CFL condition is always a few orders of magnitude greater than the other; i.e.Table 1.Parameters values.

Model implementation
We have implemented the model presented here using compute unified device architecture (CUDA), a software and hardware architecture that was developed by NVIDIA, to allow the use of their graphic cards (GPUs) as highperformance parallel-computing devices.This provides us with a speed-up of two orders of magnitude over the standard CPU-based implementation, as shown for the CUDA SPH fluid-dynamics code presented in Hérault et al. [2010], on which the present study is based.We recall here briefly the most important aspects of the algorithm implementation, referring the reader to Hérault et al. [2010] for further details.Differences with respect to this implementation are also highlighted.
The particle system is split into a CPU and a GPU part.The CPU side acts as a controller: it takes care of data initialization (defining the material properties, loading the terrain DEM onto the GPU), data retrieval during the simulation (for visualization or on-disk storage), and initiating of each time-step procedure.The GPU side holds all of the data and runs all of the computations during the simulation, with each time step decomposed into three main substeps: neighbor search, force computation, and system update.

Neighbor search and list update
The compact support of the smoothing kernels used in SPH guarantees that when computing forces, only the particles within the support provide a nonzero contribution.We can therefore limit interactions to these neighboring particles, reducing an N × N problem to a problem of the order of O(N)3 , which is more efficient if the neighbor search itself can be limited to O(N).
To achieve this, we divide the space into a regular grid where the cells have the size of the support of the kernel, kh.The cells are then linearly indexed and we sort the particles by their cell indices.A neighbor list for each particle is then built from the particles in the same cell or in adjacent cells.Although our problem is three-dimensional, lava flow tends to have a predominantly two dimensional structure, so we only consider a grid in the xy plane, with each cell indexed disregarding its z coordinate.
During this operation, we also compute S cells , the total surface of the grid cells that contain at least one particle.This surface is a good approximation of the ground contact surface of the flow, S g , which is needed in Equations ( 22) and ( 24).Due to the specificity of lava flow, a predominant twodimensional structure and a smooth free surface, we also consider that in the first approximation the extent of the free surface S f needed in Equation ( 23) is equal to S cells .
As the neighbor list construction is memory-intensive and low on the number of computations, it is the least efficient part of our GPU implementation.Therefore, rather than rebuilding the list at each iteration, this part of the code is only executed every N n time steps, with Nn being a userconfigurable parameter with a default value of N n = 10.

Force computation
To exploit the parallel nature of GPUs, force computation is computed per particle, meaning in particular that each particle-particle interaction is computed twice (once for each particle involved).
During the force computation, we calculate the density derivative and the particle acceleration from the Navier-Stokes equations in the forms of Equations ( 17) and ( 18), as described earlier, and the thermal evolution from the equations in Section 4.3.
In addition to the ground interaction described in Section 4.4, particles in the vent area are also subject to an additional repulsive force that lifts them, which ensures that new particles can be injected underneath.The lift force is computed so as to ensure that the injection of new particles can be carried out at a rate that matches the simulated effusion rate of the vent.
The force computation routine also returns the maximum allowed time step according to the CFL condition discussed in Section 4.5.For nonNewtonian rheologies, we also have to compute the apparent viscosity of each particle before any force computation.

System update
After force computation, the system is updated using the numerical integration scheme described in Section 4.5.If the lava effusion rate determines that new material is to be input into the system, new particles are generated at the vent location.
To determine the frequency with which new particles are injected, the system precomputes the volume corresponding to one layer of particles that occupy the vent area.The ratio of this volume to the flux rate determines the time at which the next injection will occur.After every injection, a neighbor list construction will be forced even if less than the predefined N n steps have occurred, to ensure that the new particles will properly interact with the rest of the system.

Results
Due to the lack of sufficiently complete and detailed information on the eruptive events on Mount Etna, no comparisons against these have been done yet.Instead, a qualitative comparison with a variety of combinations of parameters was carried out to assess the influence that the different rheologies and the temperature dependency have on the emplacement of the flow according to the model.For all of the simulations we used the parameters listed in Table 1.
Simulations were carried out for each of the modeled rheologies (Newton, modified Bingham model, and modified Herschel-Bulkley model), both in the case of temperature-independent parameters and in the case of temperature-dependent viscosity, and (for a modified Bingham model and a modified Herschel-Bulkley model), yield strength.Temperature-dependent examples are further divided into simulations with a higher solidus temperature, higher viscosity and higher flux rate, versus simulations with lower solidus temperature, lower viscosity and both higher and lower flux rates.
The temperature-dependent viscosity uses the law from Giordano and Dingwell [2003]: where T is the temperature in K, and [H 2 O] is the water content in weight percentage.Higher water contents result in flows with lower viscosity, and vice versa.
For the yield strength, we use [Miyamoto and Sasaki 1997]: Simulations were carried out using the real topography from Mount Etna, for the south summit zone.The DEM for the simulation was derived from a 2 m resolution DEM from 2005, interpolated to a 1 m resolution.
In all of the figures, the kinematic viscosity is expressed in m 2 • s −1 , the yield strength in Pa, the velocity in m • s −1 , the temperature in K, and the shear rate in s −1 .In the Figures containing surface temperatures, we plotted the temperature drop from the vent temperature.

Simulation #1: Newtonian fluid
The first simulation results that are shown in Figures 5  and 6 uses the following parameters: Newtonian fluid, constant kinematic viscosity o = 110 m 2 • s −1 , T s = 1263 K, and q = 25 m 3 • s −1 .
The discretization uses Dp = 1 m, which results in about 250,000 particles after 166 min.

Simulations #2 to #4: temperature-independent Newtonian model, modified Bingham model, and modified Herschel-Bulkley model
Figures 9, 10, 11 and 12 show comparisons between the evolution of the fluid with different rheology laws.The discretization still uses Dp = 1 m, which results in about 105,000 particles after 70 min.In all three of these cases, the same solidus temperature, T s = 1263 K, and flux rate, q = 25 m 3 • s −1 , were used.The comparison involves: • a Newtonian fluid with o = 14.5 m 2 • s −1 ; • a modified Bingham model with o 0 = = 10 m 2 • s −1 , x 0 = 300 Pa, and m = 1,000; the apparent viscosity o app = , with n app computed from 14, is show in Figure 7; • a modified Herschel-Bulkley model with o 0 = = 10 m 2 • s −1 , x 0 = 300 Pa, t 1 = 0.01 s, and m = 1,000; the apparent viscosity o app = , with n app computed using Equation ( 13), and is plotted in Figure 8.
With these parameters, the modified Herschel-Bulkley model represents a shear thinning fluid, with an apparent viscosity that decreases as the shear rate increases.The strongest differences between the two models are to be expected in these cases, in particular in the areas where the fluid velocity increases due to variations in the slope.This behavior can indeed be observed towards the final part of the simulation, as the lava gets funneled into a depression of the topography (Figure 12).Again, we compare the Newtonian model with the modified Bingham and modified Herschel-Bulkley models, with the m and t 1 parameters as in Section 6.2.The temperature-dependent viscosity parameters o(T) (Newton), x(T) and o app (T) (modified Bingham and modified Herschel-Bulkley models) are plotted in Figures 13,14,15 and 16,respectively.In Equations ( 13) and ( 14) we take n 0 = to(T) and x 0 = x(T).
Compared to the previous simulations (Figure 12), the differences between the modified Bingham and modified Herschel-Bulkley models grow more significant as the lava progresses in the funneling depression provided by the topography (Figure 20).
In the Newtonian case, the lava overflows the computational domain before t = 70 min (Figures 17 and 18).

Simulations #8 and #9: higher resolution Newton and modified Bingham models, with higher water content
The results presented here are shown in Figures 24 and 25, and they use a higher resolution, with Dp = 0.5 m, which results in about 230,000 particles after 20 min.
The dynamic viscosity now depends on T, using Equation ( 26), with [H 2 O] = 0.02; and the yield strength follows Equation ( 27).The Newtonian and modified Bingham models are compared, with T s = 1263 K, q = 25 m 3 • s −1 and                [H 2 O] = 0.02.For Bingham we have m = 500.The temperature-dependent functions o(T), x(T) and o app (T) for this case are plotted in Figures 21, 22 and 23.Again, in Equations ( 13) and ( 14), we take n 0 = to(T) and x 0 = x(T).

Simulations #10 and #11: lower effusion rates
These results are shown in Figures 26 and 27, and they where obtained with the same parameters as the previous simulation, but with a flux rate q = 10 m 3 • s −1 , which results in about 95,000 particles after 20 min.

Comments on the simulations
Qualitatively, the simulation results confirm that neither the Newtonian rheology, even with temperature dependency, nor the temperature-independent, nonNewtonian rheologies are compatible with the Etnean lava flow.The results thus confirm that a nonNewtonian model with temperaturedependent parameters must be used, as only in these cases can the flow develop levees and channels, such as those that are found on Etna, as seen in Figure 28.

Conclusions
A fully three-dimensional SPH lava-flow model is presented in this study.The classic SPH model for fluid dynamics is integrated with the thermal model of Cleary and Monaghan [1999] and a phase transition model similar to that of Monaghan et al. [2005]; a parameter is used to model the fraction of latent heat lost (or gained) during phase transition, allowing us to easily track the solidification fronts in the flow.
Following recent developments in the physical modeling of lava-flow rheology, both Bingham and power-law fluids can be simulated by our SPH code.The topography of the terrain is provided as a DEM, and we use the GPU built-in hardware interpolation to obtain the ground height and the surface normal, at any point of the surface.This allows us to implement the boundary conditions using Lennard-Joneslike repulsive forces normal to the terrain at the particle location, without actually covering the domain boundary with particles.In addition to being computationally cheaper, our approach ensures that a single plane exerts a constant force on a particle moving parallel to it, which solves one of the typical issues of repulsive particles.Our model also includes a limited representation of particle aggregation, which is restricted to the accretion of the topography from solid lava.An appropriate extension of this representation is important for the correct modeling and treatment of crust formation.
A number of example simulations are presented in the present study, with a variety of physical models and rheological parameters.Although no comparisons are carried out against real events, the results from the simulations are sufficient to strongly highlight the importance of the rheology in the flow emplacement, both in terms of the actual physical model (Newton model vs. modified Bingham model vs. modified Herschel-Bulkley model with a power different from 1) and in terms of the parameter values.
The simulation times are currently very long, and even with the benefit of the highly parallel GPU hardware, the computational load restricts the application to relatively short simulations.
To extend the applicability of the model to the simulation of complete eruptions, which would include both scenario forecasting and validation against case studies, the following improvements are necessary: • The computational run-time must be lowered.
• A detailed and accurate rheological model of lava must be obtained, including both the exact physical model of the fluid and the laws of the temperature dependency for the rheological parameter.
• Accurate and frequent effusion rate estimates must be available for the simulation of real events.
To reduce the simulation run-time, we are focusing our research on the distribution of the computation across multiple GPUs, to exploit the strong parallel nature of the SPH method.
A better assessment of the input parameters to the model (rheology model, rheological parameters, and effusion rates), on the other hand, is a rather more complex issue that also requires more accurate field measures during eruptions, as well as a better understanding of the physical model to be used for lava.For this, our model can actually be used for a reverse validation.
Indeed, the nonNewtonian models as implemented in our SPH code can be validated against other fluids with known rheology and thermal dependency.After the validation of the model, experimental data obtained by melting lava samples in controlled environments at different temperatures can be compared with a wide range of numerical experiments in an inverse problem, to determine the rheological model and parameter laws that best fit the experimental data.
A similar approach can be used for the estimation of the flux rates.An approach that has gained a growing interest uses infrared satellite imagery [Harris et al. 1998] to reconstruct the mass flux rates from the thermal flux.This approach has some well-known limitations, and with a lack of other means for accurate and frequent estimation of the flux rates during eruptive events, the only possible validations are a posteriori, by comparing the estimated total volume emitted with the final lava emplacement.
Once again, our model can be used in an inverse problem to determine the actual error caused by the application of the method from Harris et al. [1998], and can possibly be used to determine more accurate laws that correlate thermal and flux rates for lava flow.This can be achieved by running a number of simulations with different constant and variable flux rates, and determining the evolution of the surface radiance of the lava emplacement during the course of an eruption.

Figure 1 .
Figure 1.Interpolation of normal to the ground.

Figure 2 .
Figure 2. The DEM data (blue) and the terrain normals (red) generated at a resolution ten times less than that used for the DEM.

Figure 3 .
Figure 3.A particle in contact with the ground.

Figure 4 .
Figure 4. Free surface detection.(a) Free surface particle.(b) Particle that is in the internal part of the flow.
modified Bingham model, 18, 19 and 20  show the results for the temperature-dependent viscosity parameters.Again, the discretization uses Dp = 1 m, leading to about 105,000 particles at 70 min, with the solidus temperature T s = 1263 K, and a flux rate q = 25 m 3 • s −1 .The water content for the temperature-dependent viscosity and yield strength (Equations 26 and 27) is set to [H 2 O] = 0.005.

Figure 5 .
Figure 5. Evolution of temperature (left; temperature drop from vent, K) and solidification (right), at

Figure 6 .
Figure 6.Evolution of velocity field (left; m • s −1 ) and particle aggregated to the ground (right), at t = 8.3

Figure 9 .
Figure 9. Evolution of apparent viscosity (m 2• s −1 ), as a comparison between the Newtonian (left) and

Figure 18 .
Figure 18.Evolution of apparent viscosity (m 2 • s −1 ) as a comparison between the Newtonian (left) and modified Bingham (right) models with temperature dependent viscosity and yield strength, at t = 10 min (top), t = 40 min (middle) and t = 70 min (bottom).The missing 70-min panel for the Newtonian model (bottom, left) arises as the lava overflows the computational domain.

Figure 23 .Figure 24 .
Figure 24.Evolution of temperature drop from vent (K) as a comparison between the Newtonian (left) and modified Bingham (right) models with temperature-dependent viscosity and yield strength, at t = 7

Figure 25 .
Figure 25.Evolution of apparent viscosity (m 2 • s −1 ) as a comparison between the Newtonian (left) and modified Bingham (right) models with temperature-dependent viscosity and yield strength, at t = 7 min

Figure 26 .
Figure 26.Evolution of temperature drop from vent (K) as a comparison between the Newtonian (left) and modified Bingham (right) models with temperature-dependent viscosity and yield strength, at t = 7

Figure 27 .
Figure 27.Evolution of apparent viscosity (m 2 • s −1 ) as a comparison between the Newtonian (left) and modified Bingham (right) models with temperature-dependent viscosity and yield strength, at t = 7 min