2.1 Description of the model

During an explosive eruption, a volcanic jet erupts out from a vent with a speed of several tens to more than 150 m s⁻¹, driven by expanding gas. The jet is initially denser than the surrounding atmosphere and begins to decelerate through negative buoyancy and turbulent interaction with surrounding air. Cauliflower-like vortices are generated along jet margins, within which process, air is entrained and heated up, reducing the bulk density of the entire jet, in many cases, to less than that of the surrounding atmosphere. Once it becomes buoyant, such a jet develops into a Plinian or sub-Plinian plume, rising up to several kilometers to tens of kilometers until its heat is exhausted. Jets that lose their momentum before becoming buoyant collapse back onto the ground and transform into pyroclastic flows, surges and ignimbrites. During the process of plume rising up, relatively larger particles might separate from main stream of the plume, falling down onto the ground and possibly be re-entrained into the plume at a lower height (Ernst et al., 1996). Within this process, erupted vapor condenses to liquid (droplet) and even further to ice. Latent heat released from phase change of erupted vapor further heats up entrained air and further dilutes the bulk density. The entrained vapor might also experience a similar process and impact plume development. Particle aggregation processes (Carey and Sigurdsson, 1982; Taddeucci et al., 2011), either due to presence of liquid water, resulting from particle collision or driven by electrostatic forces, might occur inside plume and thereby affect the sedimentation.

All in all, the process of plume development is essentially a multiphase turbulent mixing process coupled with heat transfer and other microphysical and chemical reactions.

As an initial effort towards exploiting the feasibility and advantage of SPH in plume modeling, our model is designed to describe an injection of well-mixed solid and volcanic gas from a circular vent above a flat surface into a stratified stationary atmosphere following SK-3D (Suzuki et al., 2005). In this model, molecular viscosity and heat conduction are neglected since turbulent exchange coefficients are dominant. Erupted material consisting of solid materials with different size and mixture of gases is assumed to be well mixed and behave like a single-phase fluid (phase 2) which is valid for eruptions with fine particles and ash. Air (also a mixture of different gases) is assumed to be another phase (phase 1). Thermodynamic equilibrium is assumed so that no separate energy equation is needed for each phase. As a result, there is only one energy equation for both phases (heat exchange term between different phases does not show up under this assumption). Dynamic equilibrium is assumed so that no separate momentum equation is needed for each phase. As a result, there is only one vector momentum equation for both phases (drag force term does not show up with this assumption).

Because of the above assumptions, all other microphysical processes (such as the phase changes of H₂O, aggregation, disaggregation, absorption of gas on the surface of solids, solution of gas into a liquid) and chemical processes are not considered in this model. These ignored microphysics factors would play critical roles under particular eruption conditions. Capturing of these processes needs a more comprehensive model. One critical element in plume development, the effect of wind, is also not yet considered in this model. Introducing wind effects in our model requires dynamic pressure boundary conditions, which requires more numerical effort and algorithm design. To summarize, our model is not valid for eruptions where wind effect plays a significant role in its development, usually referred to as a weak plume. Our model also lacks the ability in modeling plumes with large particles or eruptions in which microphysics plays non-ignorable roles, such as an eruption of the El Chichón volcano on 4 April 1982 (Sigurdsson et al., 1984; Folch et al., 2016). We are focused here on developing the SPH-based methodology in the context of the more basic (and more critical) aspects of the volcanic plume and therefore devote our effort to this relative simpler model. It is worthwhile to mention here that because SPH is adopted as our numerical method, adding these physics into our model would require much less work in terms of programming compared to mesh-based methods. Since our plan is an open-source distribution of the tool, we believe some of these enhancements will rapidly ensue with community participation.

2.2 Governing equations

Based on above assumptions, the governing equations of our model are given as (which is the same as the governing equations of SK-3D) (Suzuki et al., 2005)

where ρ is the density, v is the velocity, ξ is the mass fraction of ejected material, g is the gravitational acceleration, and I is a unit tensor. $E=e+K$ is the total energy which is a summation of kinetic energy K and internal energy e. An additional equation is required to close the system. In this model, the equation for closing the system is the following equation of state (EOS):

where

where C_v is the specific heat with constant volume, and R is the gas constant. ξ is the mass fraction of erupted material. The subscript “m” represents mixture of ejected material and air, “s” represents solid portion in the ejected material, “g” represents gas portion in the ejected material, “a” represents air, and 0 represents physical properties of erupted material. ξ_g0 is the mass fraction of vapor in the erupted material.

In mesh-based methods, governing equations in Eulerian form, Eqs. (1)–(4), are directly discretized. For SPH, governing equations in Lagrangian form are needed. By deducting kinetic energy from the energy equation, subtracting mass conservation from the momentum equation, combining the transient term and advection term into the material derivative term (for any function A, material derivative is defined as $\frac{\text{d}A}{\text{d}t}=\frac{\partial A}{\partial t}+\mathit{v}\cdot \mathrm{\nabla }A$), the governing equations are put into the final form, in which the advection term does not appear explicitly.

2.3 Boundary conditions

In the current model, the initial domain is a box. The boundaries are categorized as the velocity inlet (a circular area at the center of the bottom of the box), wall boundary (box bottom) and pressure outlet (other faces of the box).

3.1 Fundamental principles

There are several procedures for discretizing governing equations (PDEs or ODEs) with SPH. We present here one of them following Monaghan (1992, 2005, 2012). The starting point of approximating a function with SPH is the translation property of the Dirac function δ(x), for an arbitrary function A(x), the following equation holds.

The Dirac function δ in Eq. (25) can be approximated by a weighting function $w\left(\mathit{x}-{\mathit{x}}^{\mathbf{\prime }},h\right)$ (or $\left.w\left({\mathit{x}}^{\mathbf{\prime }}-\mathit{x},h\right)\right)$ which tends to a Dirac function when the smoothing length h→0:

The weighting function, as an approximate form of the Dirac function, should satisfy the normalization condition:

Besides normalization, the weighting function of particle a has to be symmetric with respect to a to ensure that neighbor particles located at the same distance away from a contribute equally to the SPH summation equation; see Eq. (28).

The weighting function also needs to satisfy conditions such as positivity and compact support. In addition, the kernel function must be monotonically decreasing with the distance between particles.

There is a wide variety of possible weighting functions that can satisfy these requirements, such as spline functions (with different orders) and Gaussian functions. Generally, the accuracy increases with the order of the polynomials of the kernel function, but the computational cost also increases as the number of interactions increases. We adopt a truncated Gaussian function as the weighting function in our simulation.

where d is the number of dimensions. The derivative of the weighting function is

By replacing the δ function in Eq. (26) with the kernel function w, an arbitrary function A(x) can then be approximated by

As the weighting function is symmetric (Eq. 28) and satisfies the normalization condition (Eq. 27), odd error terms in Eq. (31) vanish, leading to a second-order approximation. However, in practice, a second order of accuracy cannot be achieved because there is no guarantee of the symmetry of particle distribution in real simulation (Price, 2012). Recall that $\text{d}{\mathit{x}}^{\mathbf{\prime }}=\frac{\text{d}m\left({\mathit{x}}^{\mathbf{\prime }}\right)}{\mathit{\rho }\left({\mathit{x}}^{\mathbf{\prime }}\right)}$; the integration equation, Eq. (31), can be approximated by summation and can lead to an approximation of the function A:

where the summation is over all the particles within the region of compact support (see Eq. 29) of the weighting function. Gradient terms may be straightforwardly calculated by taking the derivative of Eq. (32), giving

For vector quantities, the expressions are similar, simply replacing A with A in Eq. (32) and Eq. (33), giving

3.2 Artificial viscosity

In classical SPH, shock waves are handled by introducing artificial viscosity, a term that is defined based on second derivatives of velocity, to smear out discontinuities. As in the case of first-order derivatives, second-order derivatives can be estimated by differentiating a SPH interpolation twice. However, such a formulation has two disadvantages: first, it is very sensitive to irregular distribution of particles; second, the second derivative of the kernel can change sign and lead to unphysical representations (for example, viscous dissipation causes decrease of the entropy).

One of the most commonly used models of artificial viscosity (Monaghan and Gingold, 1983) is

The coefficient ν is defined as

where

The artificial viscosity term Π_ab is a Galilean invariant and vanishes for rigid rotation. It produces a repulsive force between two particles when they are approaching each other and an attractive force when they are receding from each other.

The SPH viscosity can be related to a continuum viscosity by converting the summation to integrals (Monaghan, 2005). It has been shown that shear viscosity coefficient $\mathit{\eta }=\frac{\mathit{\rho }\mathit{\alpha }hc}{\mathrm{8}}$ and bulk viscosity coefficient $\mathit{\zeta }=\frac{\mathrm{5}\mathit{\eta }}{\mathrm{3}}$ are appropriate for two-dimensional flows. $\mathit{\eta }=\frac{\mathit{\rho }\mathit{\alpha }hc}{\mathrm{10}}$ and $\mathit{\zeta }=\frac{\mathrm{5}\mathit{\eta }}{\mathrm{3}}$ are appropriate for three-dimensional flows. An extra term was added to ν considering aspects of the dissipative term in shock solutions based on Riemann solvers and led to a new formulation of artificial viscosity. We adopt this new formulation in our simulation:

where

α and β are two parameters that can be adjusted for different cases. α=1 and β=2 are recommended by Monaghan for best results. In our simulation, these two parameters are calibrated to α=0.3 and β=0.6. η is usually taken as 0.1 to prevent singularities.

3.3 Discretization of governing equations and extensibility

The basic interpolation given in Eqs. (32)–(37) provides a general way to obtain SPH expressions of governing equations. The problem is that using these expressions “as is” in general leads to quite poor gradient estimates. Various tricks can be used to conserve linear and angular momentum and thermal energy (Monaghan, 1992). Special treatments are also needed for second-order derivative terms (Monaghan, 2005). We only refer here to one of these possible discretizations of compressible Euler equations with SPH:

where a is the SPH particle index. Π is an artificial viscosity term, which is discussed in Sect. 3.2. w_ab(h) is a concise form of $w\left({\mathit{x}}_a-{\mathit{x}}_b,h\right)$, and from here on, we will use this concise form. As a Lagrangian method, particle position is also updated at every time step.

We highlight an important feature of the SPH methodology. Adding new physics and new phases into the model is trivial in terms of discretization. For example, adding a new source (or sink) into Eq. (46), adding a drag force into Eq. (47) and adding a heat exchange term into Eq. (48) lead to the new discretization form:

where the source term $\dot{\mathit{\rho }}$ can be a “sink” of erupted vapor due to its phase change. D is a drag force coefficient. κ is the heat conduction coefficient. T is the temperature. Other physics can be added easily in a similar way. Adding these new terms leads to modification of only a few lines in the source code. The drag force term should show up when dynamic disequilibrium between different phases is considered. In that case, each phase needs one set of governing equations of Navier–Stokes type. Adding a new phase into SPH code only requires adding a few new lines for the new phase besides interaction terms introduced by the new phase.

3.4 Time step

The physical quantities (velocity, density and pressure) and particle position change every time step. The Courant condition, which is in spirit similar to the Courant condition for the mesh-based methods, is used to determine the time step Δt.

where c_a is sound speed at particle a and calculated based on heat-specific ration of the mixture γ_m (see Eq. 54).

First-order Euler integration, with CFL=0.2, is used to advance in time.

3.5 Tensile instability and corrected derivatives

The classical SPH method was known to suffer from tensile instability and boundary deficiency. Tests of the standard SPH method indicate an instability in the tensile regime, while the calculations are stable in compression. A simple example of such a test calculation exhibiting the instability involves a body which is subject to a uniform initial stress, either compressive or tensile. If the initial stress is tensile, a very small velocity perturbation on a single particle can lead to particles clumping together, forming large voids and seriously corrupting density distribution. But if the initial stress is compressive, the small velocity perturbation on a single particle cannot lead to any changes in particle distribution (Swegle et al., 1995). To address these difficulties, Chen et al. (1999) proposed a corrected SPH formulation. For the 1-D case, employing a Taylor expansion for A(x) about x_a, multiplying both sides by kernel function and then doing an integration over the domain gives

Ignoring derivative terms higher than first order and writing the integral in particle approximation form leads to

Notice that the denominator in Eq. (56) is actually the summation approximation of the left side of Eq. (27). That is to say, Eqs. (56) and (32) are the same for particles far away from boundaries as the denominator in Eq. (56) becomes 1 in that case. The first-order derivative term can be obtained in a similar way:

For problems of higher dimension, the expressions for function approximation are exactly the same as Eq. (56), even though the derivation is different. The first-order derivative can be obtained by solving a system of equations explicitly or numerically (Chen et al., 1999).

3.6 Mass fraction update

Air and erupted material are represented by two different sets of SPH particles (or discretization points) in the model. Based on assumptions we made in Sect. 2, only density needs to be updated respectively for each phase. The updating of density is exactly the same as Eq. (46) in spirit. Particles of phase 1 are not counted while evaluating density of phase 2, and vice versa. Updating of density is then based on the following discretized equations.

where the subscripts a and b represent air particles (phase 1), while i and j represent particles of erupted material. β and α represent either erupted material particles or air particles. ${\mathit{\rho }}_a^{\text{a}}$ is density of phase 1 (air). ${\mathit{\rho }}_i^{\text{sg}}$ is density of phase 2 (erupted material). $\mathit{\rho }={\mathit{\rho }}^{\text{a}}+{\mathit{\rho }}^{\text{sg}}$ is density of mixture of air and erupted material. By definition, the mass fraction is updated according to Eq. (60).

In areas far away from the interface, updating of density is exactly the same as that for single-phase flow. For example, on the right side and left side (or blue areas) in Fig. 1, where there are only air particles, Eq. (59) evaluates to zero, and total density is

which is a special case of Eq. (56). For these areas occupied by only particles of phase 2 and far away from the interface, similarly, the equation for the density update becomes

That is to say, the same density updating equation can be applied for both phases and no additional numerical treatment is needed to locate where the interface is.

Figure 1In panel (a), the blue particles (phase 1) represent air particles; the red ones (phase 2) represent erupted material. Panel (b) shows corresponding mass fraction. Mass fractions are evaluated based on Eqs. (58) and (59) without any other interface track or capture method.

Download

Interface construction will become necessary and important when attempting to include the effects of mixing by resolving the detailed interface structure and dynamics of turbulence. As a Lagrangian method, interface tracking in SPH is explicit through capturing of the locations of the particles, much simpler than Eulerian methods. The existence of complex evolving interfaces between phases presents severe challenges to conventional Eulerian grid-based numerical methods. Either interface tracking (Lagrangian) (Harlow and Welch, 1965; Wrobel and Brebbia, 1991; Cheng and Armfield, 1995) or interface capturing (Eulerian) (Hirt and Nichols, 1981; Youngs, 1982; Gerlach et al., 2006; Gopala and van Wachem, 2008) methods are used to reconstruct the flow interface of free boundary flow. High computational cost, a tendency to form numerical instabilities and the inability to track complex topological changes are the significant drawbacks of tracking techniques (Hirt and Nichols, 1981; Unverdi and Tryggvason, 1992; Anderson et al., 1998). For the interface capturing (Eulerian) method, the surrendering of surface detail before the phase transport calculation means that interface reconstruction is required between time steps to recover the interface information, which needs additional numerical effort (Hirt and Nichols, 1981; Youngs, 1982). Since SPH is able to adaptively adjust the discretization and automatically construct the interface, SPH requires less additional numerical effort for interface construction and therefore is more suitable for volcanic plume simulation.

3.7 Turbulence modeling with SPH

For high-speed shearing flow, the momentum exchange and heat transfer are dominated by turbulent fluctuations as turbulent exchange coefficients are several magnitudes larger than corresponding physical coefficients (molecular viscosity and heat conduction coefficient). In addition to momentum and energy exchange, mixing between plume and air is important in plume modeling. Quantifying these mixing processes in real implementation is challenging because of the scale disparity between the large-scale fluid motion and the diffusion processes on interface that ultimately lead to mixing. Ideally, one would like to be able to include the effects of mixing on the large-scale dynamics without resolving the detailed interface structure and dynamics of turbulence to reduce computational cost. To resolve all turbulent exchange at all different scales and sub-particle-scale mixing with relative coarse resolution, a SPH sub-particle-scale (SPH-SPS) turbulence models should be included. Among existing SPH-SPS turbulence models (Holm, 1999; Monaghan, 2002; Violeau and Issa, 2007; Monaghan, 2011), we adopt a LANS-type turbulence model, the SPH−ε turbulence model (Monaghan, 2011). However, the SPH−ε turbulence model was proposed only for incompressible flow. In the following section, we will extend it for compressible flow. It is necessary to mention that all other existing SPH-SPS turbulence models (Holm, 1999; Monaghan, 2002; Violeau and Issa, 2007) also only focus on incompressible flow.

3.7.2 Turbulent heat transfer

We adopt the Reynolds analogy to get the heat transfer coefficient due to turbulence. The Prandtl number is defined as

$$\begin{array}{}\text{(73)}& \mathit{Pr}={\displaystyle \frac{C_{\text{p}}\mathit{\mu }}{\mathit{\kappa }}},\end{array}$$

where μ is the dynamic viscosity, and κ is the thermal conductivity. In addition, μ can be written in terms of the absolute viscosity (kinematic viscosity) as

$$\begin{array}{}\text{(74)}& {\displaystyle }\mathit{\mu }=\mathit{\rho }\mathit{\nu }.\end{array}$$

Then,

$$\begin{array}{}\text{(75)}& \mathit{\kappa }={\displaystyle \frac{C_{\text{p}}\mathit{\mu }}{\mathit{Pr}}}.\end{array}$$

The typical value of Pr_t for air is 0.7–0.9. We take Pr_t=0.85 for gases as recommended Kays (1994) from summarizing experimental results.

Monaghan (2005) summarized the simulation of viscosity and heat conduction in his review on SPH. We will refer to his summary in our following discussion. The additional term in the discretized momentum equation, Eq. (71), is the turbulent shear stress term. Recall that molecular viscosity can be discretized with SPH as shown in Eq. (38). It has been shown that the discretized molecular viscosity has both bulk viscosity and shear viscosity, where shear viscosity coefficient is (Monaghan, 2005)

$$\begin{array}{}\text{(76)}& {\mathit{\nu }}_t=S\mathit{\nu },\end{array}$$

with

$$\begin{array}{}\text{(77)}& S=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{10}}& \text{if}d=\mathrm{3}\\ \\ \frac{\mathrm{1}}{\mathrm{8}}& \text{if}d=\mathrm{2}\end{array}\right..\end{array}$$

The turbulent viscosity coefficient can be inferred from that formulation if we can reformulate the turbulent shear stress term in a form which is similar to the molecular shear term. Reformulating the turbulent shear stress term,

$$\begin{array}{ll}\text{(78)}& {\displaystyle }& {\displaystyle }\sum _b{\displaystyle \frac{\mathit{\epsilon }}{\mathrm{2}}}{\displaystyle \frac{m_b}{{\mathit{\rho }}_b}}{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}{\mathrm{\nabla }}_a{G}_{ab}\left(l_a\right)={\displaystyle }& {\displaystyle }\sum _b{\displaystyle \frac{\mathit{\epsilon }}{\mathrm{2}S}}{m}_b{\displaystyle \frac{{\mathit{v}}_{ab}}{{\mathit{\rho }}_b}}{\displaystyle \frac{S{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}{x_{ab}^{\mathrm{2}}}}{\displaystyle \frac{x_{ab}^{\mathrm{2}}}{{\mathit{x}}_{ab}}}{\mathrm{\nabla }}_a{G}_{ab}\left(l_a\right),\end{array}$$

the turbulent viscosity coefficient can be inferred from Eq. (78).

$$\begin{array}{}\text{(79)}& {\mathit{\nu }}_t={\displaystyle \frac{\mathit{\epsilon }}{\mathrm{2}S}}{\displaystyle \frac{{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}{\mathrm{1}}}\end{array}$$

Please note that the turbulent viscosity term has the opposite sign of the molecular viscosity term in the discretized momentum equation and there is a minus sign in the expression of Π_ab, and they cancel out.

However, the above equation is correct only for the 1-D situations. For 2-D or 3-D, it is not easy to get an explicit expression. We adopt an alternative way: obtaining a value for each pair of particles instead of persisting on getting an analytical expression. Choosing the smoothing function to be the same as the SPH kernel and the smoothing length scale l to be the same as the smoothing length h, the ratio between turbulent shear stress and physical shear stress is

$$\begin{array}{ll}\text{(80)}& {\displaystyle }{\mathrm{{\rm Y}}}_{ab}& {\displaystyle }={\displaystyle \frac{\frac{\mathit{\epsilon }}{\mathrm{2}}\frac{{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{{\mathit{\rho }}_b}}{\frac{S\mathit{\nu }}{{\mathit{\rho }}_{ab}}\frac{{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}{x_{ab}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{2}}{h}_{ab}^{\mathrm{2}}}}}{\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathit{\epsilon }\left(x_{ab}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{2}}{h}_{ab}^{\mathrm{2}}\right)}{\mathrm{2}S\mathit{\nu }}}{\displaystyle \frac{{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}.\end{array}$$

Υ_ab is essentially equivalent to the ratio between the turbulent viscous effect of particle b on particle a and molecular viscous effect of particle b on particle a. Turbulent viscosity can be easily obtained by

$$\begin{array}{ll}\text{(81)}& {\displaystyle }{\mathit{\nu }}_{t,ab}& {\displaystyle }=\mathit{\nu }{\mathrm{{\rm Y}}}_{ab}{\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathit{\epsilon }\left(x_{ab}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{2}}{h}_{ab}^{\mathrm{2}}\right)}{\mathrm{2}S}}{\displaystyle \frac{{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}.\end{array}$$

The corresponding turbulent thermal conductivity should be

$$\begin{array}{}\text{(82)}& {\mathit{\kappa }}_{t,ab}={\displaystyle \frac{\mathit{\epsilon }\stackrel{\mathrm{‾}}{C_{\text{p},ab}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}\left(x_{ab}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{2}}{h}_{ab}^{\mathrm{2}}\right){\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{\mathrm{2}S{\mathit{Pr}}_t{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}.\end{array}$$

$\stackrel{\mathrm{‾}}{C_{\text{p},ab}}$ and $\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}$ are simply the arithmetic means of specific heat and density. The term used to prevent singularity now can be removed.

$$\begin{array}{}\text{(83)}& {\mathit{\kappa }}_{t,ab}={\displaystyle \frac{\mathit{\epsilon }\stackrel{\mathrm{‾}}{C_{\text{p},ab}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}{x}_{ab}^{\mathrm{2}}{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{\mathrm{2}S{\mathit{Pr}}_t{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}\end{array}$$

We also need to prevent singularity, so

$$\begin{array}{}\text{(84)}& {\mathit{\kappa }}_{t,ab}=\left\{\begin{array}{ll}\mathrm{0}& \text{if}{\mathit{v}}_{ab}=\mathrm{0}\text{or}{\mathit{x}}_{ab}=\mathrm{0}\\ \frac{\mathit{\epsilon }\stackrel{\mathrm{‾}}{C_{\text{p},ab}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}{x}_{ab}^{\mathrm{2}}{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{\mathrm{2}S{\mathit{Pr}}_t{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}& \text{otherwise}\end{array}\right..\end{array}$$

Figure 2A cross-section view of the simulation domain in the y−z plane at 66 s. Panel (a) shows the mass fraction. Panel (b) shows all boundary conditions: the dark blue region is occupied by eruption ghost particles with a “ghost particle ID” of 0; the light blue area is occupied by pressure ghost particles with a “ghost particle ID” of 1; the gray area is filled with wall ghost particles with a “ghost particle ID” of 2. The “ghost particle ID” of all real particles is set to 100; they occupy the major portion of the domain in panel (b). Panel (c) shows the cross-section view of domain decomposition based on SFC. The simulation is conducted on 12 processors, so there are 12 subdomains in total. The cross-section view shows a portion of them.

Download

The heat conduction equation without the source term is

$$\begin{array}{}\text{(85)}& C_{\text{p}}{\displaystyle \frac{\text{d}T}{\text{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }\left(\mathit{\kappa }\mathrm{\nabla }T\right).\end{array}$$

The second spatial derivative can be approximated with SPH by following Monaghan (2005):

$$\begin{array}{}\text{(86)}& C_{\text{p}}{\displaystyle \frac{\text{d}T}{\text{d}t}}=\sum _b{\displaystyle \frac{m_b}{{\mathit{\rho }}_a{\mathit{\rho }}_b}}\left({\mathit{\kappa }}_a+{\mathit{\kappa }}_b\right)\left(T_a-T_b\right)F_{ab}\left(h\right),\end{array}$$

where F_ab(h) is short for $F\left({\mathit{x}}_a-{\mathit{x}}_b,h\right)$, whose definition is

$$\begin{array}{}\text{(87)}& F_{ab}\left(h\right){\mathit{x}}_{ab}={\mathrm{\nabla }}_a{w}_{ab}\left(h\right).\end{array}$$

F_ab is always non-positive, which guarantees that heat flux flows from hot to cold. Plug the turbulent thermal conductivity into the heat conduction equation:

$$\begin{array}{ll}\text{(88)}& {\displaystyle }{C}_{\text{p}}{\displaystyle \frac{\text{d}T}{\text{d}t}}& {\displaystyle }=\sum _b{\displaystyle \frac{m_b}{{\mathit{\rho }}_a{\mathit{\rho }}_b}}\left({\mathit{\kappa }}_a+{\mathit{\kappa }}_b\right)\left(T_a-T_b\right)F_{ab}\left(h\right){\displaystyle }& {\displaystyle }=\mathrm{2}\sum _b{\displaystyle \frac{m_b}{{\mathit{\rho }}_a{\mathit{\rho }}_b}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{C_{\text{p},ab}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}\mathit{\epsilon }{x}_{ab}^{\mathrm{2}}{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{\mathrm{2}{\mathit{Pr}}_{t}S{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}\left(T_a-T_b\right)F_{ab}\left(h\right).\end{array}$$

Notice that the number “2” in the front of Eq. (88) comes from integration approximation of the second-order derivative (Cleary and Monaghan, 1999). By further simplification, we get

$$\begin{array}{ll}\text{(89)}& {\displaystyle }& {\displaystyle }{C}_{\text{p}}{\displaystyle \frac{\text{d}T}{\text{d}t}}={\displaystyle }& {\displaystyle \frac{\mathit{\epsilon }}{S{\mathit{Pr}}_t}}\sum _b{\displaystyle \frac{m_b}{{\mathit{\rho }}_a{\mathit{\rho }}_b}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{C_{\text{p},ab}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{ab}}{x}_{ab}^{\mathrm{2}}{\mathit{v}}_{ab}\cdot {\mathit{v}}_{ab}}{{\mathit{v}}_{ab}\cdot {\mathit{x}}_{ab}}}\left(T_a-T_b\right)F_{ab}\left(h\right).\end{array}$$

3.8 Boundary conditions

All boundary conditions are imposed by ghost particles. Figure 2b shows how boundaries are deployed.

3.9 Parallelism and performance

One disadvantage of the 3-D model is that it usually takes a much longer time than 1-D models to complete one simulation. This disadvantage further prevents simulation with finer resolution and accounting for more physics in one model. Non-intrusive uncertainty analysis, which is commonly adopted in hazard forecasting, requires finishing multiple simulations within a given time window. High-performance computing is therefore essential. Among existing CPU parallel SPH schemes, most of them focus on the neighbor search algorithm and dynamic load balancing (Ferrari et al., 2009; Crespo et al., 2015). Less attention has been paid to developing more flexible data management schemes for more complicated problems. Motivated by techniques developed for mesh-based methods, we develop a complete framework for parallelizing SPH with distributed memory parallelism (MPI) allowing flexible and efficient data access.

The time complexity for an SPH method is O(N²), where N is the total number of SPH particles. Efficient neighbor searches and compact-supported kernel functions can help to reduce computational cost. Of the many possible choices, we adopt a background grid which was proposed by Monaghan and Lattanzio (1985) and is quite popular in parallel SPH. Then, the time complexity is reduced to O(MN+mN), where m is the average of number of particles within the compact support of the kernel function, M is the average number of particles in subdomain within which neighbor searches are carried out. This background grid is also used for domain decomposition. We refer to the elements of the background grid, namely squares for two dimensions and cubes for three dimensions, as buckets. As for the actual storage of the physical quantities, different strategies have been adopted in existing implementations of SPH.

In DualSPHysics (Crespo et al., 2015), the physical quantities of each particle (position, velocity, density, etc.) are stored in arrays. The particles (and the arrays with particle data) are reordered following the order of the cells. This has two advantages: (1) the access pattern is more regular and more efficient; (2) it is easy to identify the particles that belong to a cell by using a range since the first particle of each cell is known. But adding, deleting and especially accessing particles are cumbersome. Ferrari et al. (2009) adopted linked lists using pointers so that particles can be deleted or added during the simulation. Storage problems caused by fixed-size arrays are thereby also eliminated. We define C++ classes which contain all data of particles and buckets. As for the management of data, we adopt hash tables to store pointers to particles and buckets, which gives us not only flexibility of deleting and adding elements but also quicker access compared with linked lists. Instead of using the “natural order” to number particles, we adopt a SFC-based index to give each particle and background bucket a unique identifier – a strategy known to preserve data locality at minimal cost. The SFC-based numbering strategy is further extended to include time step information so that particles added at the same position but different time have different identifiers.

The parallelization is achieved by splitting computational domain into subdomains. Each subdomain is computed by a single processor. For any subdomain, information from its neighboring subdomains is required when updating physical quantities. To guarantee consistency, data are synchronized if physical quantities are updated. Even though more complicated graph-based partitioning tools (Biswas and Oliker, 1999) might get higher quality decomposition, they require much more effort in programming and computation. So we adopt an easy programming scheme based on SFC (Patra and Kim, 1999) to decompose the computational domain. Figure 2c shows a cross-section view of domain decomposition.

More details about the data structure, domain decomposition, load balancing, domain adjusting and performance benchmarking have been published separately (Cao et al., 2017).

Figure 3Strong scalability test result.

Download

Figure 4Weak scalability test result.

Download

Figure 5The effect of domain adjusting on simulation time. Panel (a) shows execution time without domain adjusting; panel (b) shows execution time with domain adjusting. Different bins represent execution time up to specific physical time indicated by horizontal axis.

Download

Performance tests have been carried out on the computational cluster of Center for Computational Research (CCR) at the University at Buffalo. Intel Xeon E5645 CPUs running at 2.40 GHz clock rate with 4 GB memory per core on a Q-Logic Infiniband are used in these tests. Each node is comprised of two sockets with six of these cores. Memory and level 3 cache are shared on each node. The initial domain is [−4.8 km, 4.8 km] × [−4.8 km, 4.8 km] × [0 km, 6 km]. Almost linear speed up is observed in our strong scalability test (Fig. 3).

The weak scalability test is conducted with the same initial domain and various smoothing lengths. Each simulation runs for 400 time steps. The average number of real particles of each process keeps constant at 25 900. As shown in Fig. 4, simulation times increase around one-third when the number of cores increases from 16 to 4208. For the test problem in this section, the volcanic plume finally reaches a region of [−30 km, 30 km] × [−30 km, 30 km] × [1.5 km, 40 km] after around 400 s of eruption. When numerical simulation goes up to 90 s, the plume is still within a region of [−10 km, 10 km] × [−10 km, 10 km] × [0 km, 25 km]. This implies that adjusting of domain can avoid computing a large number of uninfluenced air particles, especially for the beginning stage of simulation. A domain-adjusting algorithm (Cao et al., 2017) is designed and implemented in our code. Figure 5 shows that simulation time of the test problem is greatly reduced when we adopt the domain-adjusting strategy.

4.1 1-D shock tube tests

1-D shock tube tests are first conducted to verify our code. Input parameter of each tests can be found in Table 1. These tests can represent typical cases in 1-D. Test 1 consists of a left rarefaction, a right traveling contact and a right shock. Density decreases downwind of the contact wave. Test 2 also consists of a left rarefaction, a right traveling contact and a right shock. Density increases downwind of the contact wave. Test 3 is a double expansion test with different initial density.

Table 1Input parameters of 1-D shock tube tests.

Download Print Version | Download XLSX

In Table 1, subscript “L” refers to the left side and “R” to the right side; t_f is the total simulation time. The initial interval between two adjacent particles is 0.03. The computational domain is $[-\mathrm{0.4},\mathrm{0.4}]$. The specific internal energy is compared against exact solutions. As shown in Fig. 6, the position and magnitude of the waves are correctly predicted. The fluctuations near the contact discontinuity are caused by sharp change of smoothing length.

Figure 6Comparison of specific internal energy of simulation results against analytical results for shock tube tests. The plots from left to right correspond to test 1, test 2 and test 3, respectively.

Download

4.2 Simulation of JPUE

JPUE can be considered as a simplified volcanic plume. While the effect of stratified atmosphere and the effect of expansion due to high temperature in volcanic plume are not represented, JPUE reproduces the entrainment due to turbulent mixing which is one of the key elements in the volcanic plume development. There exist consistently good experimental data (List, 1982; Dimotakis et al., 1983; Papanicolaou and List, 1988; Ezzamel et al., 2015) that describe the JPUE flow field giving insight into details of JPUE, such as transverse velocity and concentration profile. In this section, we verify that our code and the SPH−ε turbulence model is able to reproduce the features of turbulent entrainment by a JPUE simulation.

As many of these experiments were conducted with liquid, we replace the original equation of state (Eq. 5) with a weakly compressible Tait equation of state (Becker and Teschner, 2007) (see Eq. 94) to avoid solving the Poisson equation:

with γ=7 and B is evaluated by

where c is the speed of sound in the liquid. The energy equation is actually decoupled from the momentum conservation equation and the mass conservation equation by using this EOS. In addition, the “atmosphere” is assumed to be uniform and gravity is set to be zero. We set the temperature and density of ejected material the same as surrounding ambient. This further simplifies the scenario for the convenience of studying turbulent mixing.

One overall feature of JPUE is “self-similarity”, which means that the evolution of the JPUE is determined solely by the local scale of length and velocity, which theoretically account for the fact that the rate of entrainment at the edge of JPUE is proportional to a characteristic velocity at each height. As a result, physical and numerical experiments do not necessarily have exactly the same setups and are compared on a non-dimensional basis.

Table 2List of eruption condition for the test cases.

Download Print Version | Download XLSX

Figure 7Dimensionless concentration and velocity distribution across the cross section.

Download

Figure 8Panel (a) shows normalized jet width (which determined based on velocity) along the centerline. Panel (b) shows normalized concentration along the centerline.

Download

A three-dimensional axisymmetric JPUE which ejects from a round vent is simulated with eruption parameters listed in Table 2. Material properties of water are used as material properties for both phases. The results are compared with experiments (George et al., 1977; Papanicolaou and List, 1988) for verification purposes. Experimental data of concentration and velocity distribution across the cross section were fit into a Gaussian profile (see Eq. 96) by Papanicolaou and List (1988) and George et al. (1977) even though the actual profiles are slightly different from the Gaussian profile.

where φ is either velocity or concentration, the subscript “c” represents the centerline. r is the distance from the centerline on any cross section. z is the axial distance from the origin of the jet transverse section under consideration. The coefficient “coef” for concentration is 80 and 50, respectively, according to George et al. (1977) and Papanicolaou and List (1988). The “coef” for velocity is 90 and 55, respectively, according to George et al. (1977) and Papanicolaou and List (1988).

Papanicolaou and List (1988) also fit concentration distribution and jet width based on velocity along the centerline into a straight line (see Eq. 97).

where subscript 0 represents the cross-sectionally averaged exit value; D is the diameter of vent. The “slope” for jet width based on velocity is 0.104 and for concentration is 0.157. The “intercept” for jet width based on velocity is 2.58, while that for the concentration is 4.35.

Although both velocity and concentration are found to be well matched with experimental results, a small disparity in both velocity and concentration is observed near the boundary of the jet, which is possibly caused by overestimation of the drag effect by standard SPH (Ritchie and Thomas, 2001). Ritchie and Thomas (2001) also proposed an alternative way to update density which relieved the overestimating of the drag effect. However, how well this method conserves mass is not clear. There are several other factors that might also attribute to such disparity. The Reynolds number is not reported in many experiments assuming a high enough Reynolds number. In addition, some details of the experiments, such as exit velocity profile and viscosity of the experimental liquid, are not reported. These factors prevent us from numerically reproducing these experiments in an exact way as they were conducted. However, the features of JPUE are correctly reproduced with our code.

We also investigated the mixing due to turbulence in the JPUE simulation by checking the mixture of the two phases. It is shown in Fig. 9 that the ejected material and ambient fluids are mixed through eddies at the outer shear region. Also, the inner dense core dispersed gradually due to erosion of the outer shear region. Hence, our confidence in the numerical correctness of our code is greatly reinforced.

Figure 9Panel (a) shows particle distribution. Particles of phase 1 (blue) are gradually entrained and mixed with erupted particles (red) as jet flows downstream. Panel (b) shows the mass fraction of erupted material at the moment corresponding to panel (a).

Download

4.3 Simulation of a volcanic plume

The development of a volcanic plume is more complicated than JPUE in several aspects. Besides turbulent entrainment of ambient fluids, development of the volcanic plume also involves heating of entrained air and expansion in a stratified atmosphere. A strong eruption column without wind is tested in this section for the purpose of further verification and validation. Both global variables and local variables are compared with existing models.

4.3.1 Input parameters

Eruption parameters, material properties and atmosphere are chosen to be the same as the strong plume no-wind case in a comparison study on eruptive column models by Costa et al. (2016). Eruption conditions are listed in Table 2. As our model does not distinguish solid particles of different sizes, only mass fraction of solids of all sizes is used in simulation even though two size classes were provided by Costa et al. (2016). The density of erupted material at the vent and radius of the vent can be computed from the given parameters. The eruption pressure is assumed to be the same as the ambient pressure at the vent and hence is not given in the table. The vertical profiles of atmospheric properties were obtained based on the reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) for the period corresponding to the climactic phase of the Mt. Pinatubo eruption (Philippines, 15 June 1991). These conditions are more typical of a tropical atmosphere (see Fig. 1b in Costa et al., 2016). Vertical distribution of temperature, pressure and density is used to establish stratified atmosphere. Wind velocity and specific humidity are not used in our simulation even though they were also provided by Costa et al. (2016) (see Fig. 1b). Material properties, shown in Table 3, are selected based on properties of the Pinatubo and Shinmoedake eruptions. Other material properties not given in the table can be inferred from these given parameters based on their relationships.

Table 3List of material properties.

Download Print Version | Download XLSX

Figure 10Mass fraction for t=500 s after eruption. Panels (a, b) are visualizations of SPH simulation results. Panel (b) shows visualization of a slice of the computational simulation, whose thickness is around 10 000 m. The lowest portion of the plume represents erupted material in the eruption vent (the underground portion). Panels (c, d) show the contour of the mass fraction and velocity quiver on an x−z plane. These panels are plotted utilizing postprocessed data (see Appendix A). The contour levels in the plot are 0.00001, 0.0001, 0.001, 0.1, 0.3, 0.75 and 0.95. Panel (d) is a zoomed view of velocity quiver showing plume expansion and entrainment of air.

Download

Figure 10 shows the mass fraction of the simulated volcanic plume at 500 s after eruption, at which time the plume starts spreading radially. A contour plot of the mass fraction on the vertical cross section (x−z plane) was also provided. The zoomed view of the quiver plot shows detailed entrainment of air at the margin of the plume.

4.3.2 Global and local variables

One of the key global quantities of great interest is the altitude to which the plume rises. The top height predicted by our model is around 40 km which agrees with other plume models. For example, the height predicted by PDAC is 42 500 m, by SK-3D is 39 920 m, by ATHAM is 33 392 m and by ASHEE is 36 700 m. As for local variables, the profiles of integrated temperature, density, mass fraction of entrained air, gas mass fraction, mass fraction of solid materials and the radius of the plume as a function of height are compared with existing 3-D models in Figs. 11–14. To get rid of significant fluctuations in time and space, we conducted a time-averaging and spatial integration of the dynamic 3-D flow fields by following Cerminara et al. (2016b).

Figure 11Temperature as a function of height.

Download

Figure 12The mass fraction of entrained air, gas and solids as a function of height.

Download

As particles distribute irregularly in the space in SPH simulation results, we need to project simulation results (on irregular particles) onto a predefined grid before doing time-average and spatial integration. See Appendix A for more details of postprocessing.

The profiles of local variables match well with simulation results of existing 3-D models in a general sense. The basic phenomena in volcanic plume development are correctly captured by our model.

Figure 13Density of the strong plume without wind after reaching its top height.

Download

As the height increases, the amount of entrained air also increases. Around the neutral height, where the umbrella expands, the entrainment of air shows a slight decrease due to lack of air surrounding the column at that height. The profile for gas, which accounts for both air and vapor, shows a very similar tendency to that of entrained air. Recall that vapor condensation is not considered in our model. In addition, we assume that erupted material behaves like a single-phase fluid. So the mass fraction of gas is simply a function of entrained air (Eq. 98).

$$\begin{array}{}\text{(98)}& {\mathit{\xi }}_a+{\mathit{\xi }}_{\text{g}}={\mathit{\xi }}_a+\left(\mathrm{1}-{\mathit{\xi }}_a\right){\mathit{\xi }}_{\text{g0}}\end{array}$$

Among these 3-D models, ATHAM takes vapor condensation into account and Eq. (98) does not hold for ATHAM. However, the profile of entrained air and profile of gas predicted by ATHAM are still very close to each other, which implies that ignoring water phase change is a valid assumption for eruptions similar to this test case (strong plume with erupted water fraction in erupted material less than 5 %). This observation can be explained by the fact that air occupies a larger portion of the gas, and ignoring phase change of vapor (which is only a small portion of gas) has a slight influence on plume development. As for mass fraction of solids, similarly, Eqs. (99) and (100) hold for our model.

$$\begin{array}{}\text{(99)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_{\text{s}}=\left(\mathrm{1}-{\mathit{\xi }}_a\right)\left(\mathrm{1}-{\mathit{\xi }}_{\text{g0}}\right)\text{(100)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_{\text{s}}=\mathrm{1}-\left({\mathit{\xi }}_a+{\mathit{\xi }}_{\text{g}}\right)\end{array}$$

Figure 14Radius of the strong plume without wind after reaching its top height.

Download

PDAC, which treats particles of two different sizes as two separate phases, predicted a similar mass fraction profile. That implies that assumption of dynamic equilibrium in our model is at least valid for eruptions similar to the test case.

With more cool air entrained into the plume and mixing with the hot erupted material, the temperature of the plume decreases as the height increases as shown in Fig. 11. Meanwhile, bulk density decreases due to entrainment and expansion (Fig. 13).

Our model adopts the same assumptions and governing equations as SK-3D. However, there is still an obvious disparity between the profiles of local variables of our model and SK-3D. One of the big differences between these two models is that we adopt a LANS type of turbulence model while SK-3D adopts a large eddy simulation (LES) turbulence model. This implies that choice of turbulence model might play a critical role in plume simulation.