Control of Vent Geometry on the Fluid Dynamics of Volcanic Plumes: Insights From Numerical Simulations

We present three‐dimensional numerical simulations of eruption clouds from circular to linear fissure vents to investigate the control of vent shape on the height and stability of volcanic plumes during large explosive eruptions. Our results show that clouds ejected from circular or low‐aspect‐ratio (nearly square‐like) fissure vents can be associated with radially suspended flow (RSF) at the top of the jet region, whereas those emitted from narrow‐fissure vents are not. Non‐RSF plumes are more stable than those associated with RSF because the highly concentrated parts of the ejected mixture are easily dissipated and mixed with air near the vent. Plume height in the RSF regime decreases while that in the non‐RSF regime increases with increasing aspect ratio, even for a fixed magma flow rate. These observations suggest that the efficiency of air entrainment is influenced by the vent shape, which in turn controls the dynamics of eruption plumes.


Introduction
Large explosive volcanic eruptions can produce Plinian-type buoyant plumes and/or pyroclastic flows. The hot mixture of volcanic gas and solid pyroclasts that is ejected from the vent entrains ambient air owing to turbulent mixing, significantly reducing its density because the mixed air is heated and expands. When the erupted mixture entrains a small amount of air, the eruption cloud (i.e., the mixture of the fragmented magma, lithics, gas, and entrained air) remains heavier than the ambient air and therefore collapses, resulting in pyroclastic flows. In contrast, when the erupted mixture entrains a large amount of ambient air, the eruption cloud becomes lighter than the surrounding atmosphere and rises, thereby developing a buoyant plume. In the latter scenario, the height of Plinian-type plumes is controlled by the balance between its thermal flux at the vent and the work for lifting up the eruption mixture (and the air from lower altitudes). Therefore, whether the eruption produces buoyant plumes and/or pyroclastic flows, as well as the maximum height of the plume, are strongly controlled by the amount of air entrained by turbulent mixing (e.g., Cerminara, Esposti Ongaro, & Neri, 2016;Jessop & Jellinek, 2014;Jessop et al., 2016;Trolese et al., 2019;Woods, 1988).
Some previous studies (Glaze et al., 2011(Glaze et al., , 2017Stothers, 1989;Woods, 1993) investigated the effects of vent shape (i.e., circular and linear fissure vents) using simple one-dimensional (1D) models for buoyant plumes © 2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2020GL087038
Key Points: • Three-dimensional numerical simulations reveal the effects of volcanic vent shape on plume stability and plume height • Plumes from narrow fissures are more stable than those from circular vents owing to more efficient dissipation mechanisms near the vent • Plume height has complex dependency on vent aspect ratio because of changes in the efficiency of air entrainment Supporting Information: • Supporting Information S1 • Figure S1 • Figure S2 • Figure S3 • Figure S4 Correspondence to: Y. J. Suzuki, yujiro@eri.u-tokyo.ac (e.g., Morton et al., 1956;Woods, 1988), suggesting that the critical conditions to sustain buoyant plumes is controlled by fissure width. Such studies assume that the amount of entrained air is proportional to characteristic plume velocity at each height and that the efficiency of entrainment (i.e., the entrainment coefficient) has the same constant value for plumes from both circular and fissure vents. However, the entrainment coefficient can be affected by the different turbulent features originated from the different shape of source (Deo et al., 2007a(Deo et al., , 2007bJessop et al., 2016;Mi et al., 2005;Quinn, 1992;Zaman, 1999) and the local buoyancy at each plume height (Carazzo et al., 2006;Kaminski et al., 2005;Paillat & Kaminski, 2014). In the case of volcanic plumes, the flow features and the local buoyancy are expected to depend on each other. The flow structures, their turbulent features, and the resultant plume heights are controlled by the reversing buoyancy (i.e., from negative to positive buoyancy; e.g., Suzuki, Costa, & Koyaguchi, 2016), whereas the magnitude of buoyancy is determined by three-dimensional (3D) structures of volcanic plume. In particular, when the vent radius is large, the dense unmixed part of eruption cloud can generate "fountain-like structure," known as radially suspended flow (RSF; Neri & Dobran, 1994;Suzuki & Koyaguchi, 2012;Suzuki, Costa, Cerminara, et al., 2016), and its flow structure largely changes stability of eruption column (Suzuki & Koyaguchi, 2012). In order to understand the effects of vent shape on plume dynamics, therefore, it is required to reproduce both of the reversing buoyancy and the 3D flow structures of volcanic plumes.
Examples of large explosive eruptions from linear fissure vents vary from the 1886 eruption from a 17-km-long fissure in Tarawera, New Zealand (Sable et al., 2006;Walker et al., 1984) to the~100-km-long fissure estimated for the Youngest Toba Tuff supereruption (Costa & Marti, 2016). In this study, we extend our previous studies (Costa et al., 2018;Suzuki & Koyaguchi, 2012) that investigated the dynamics of plumes from circular vent for large eruptions and reproduce the fluid dynamics of large eruption clouds from circular to fissure vents using 3D numerical models. We particularly focus on the effects of vent geometry on plume height and on the conditions for column collapse. We show that the amount of entrained air depends on vent shape, which affects both the peripheral surface of the plume that is first in contact with the ambient air (i.e., the "entrainment area") and the effective value of the entrainment coefficient (i.e., the "entrainment efficiency") associated with different flow structures. The combination of these two effects have a fundamental control on the eruption column dynamics.

Methods and Model Parameters
To investigate the effects of vent shape on eruption cloud dynamics, we numerically simulated circular and fissure vents with different aspect ratios. To quantify the effects of vent shape, we introduced two dimensionless parameters. The first is the aspect ratio of the maximum and minimum dimensions of the vent, AR = L max /L min , where L max and L min are the fissure length and width, respectively (following the definition in Mi et al., 2005). The second reflects the deviation of the vent shape from a true circle, γ = 4πA/P 2 , where A and P are the vent area and perimeter, respectively. For a given shape (circular or rectangular), γ is a function of AR; γ and AR have a maximum and minimum value, respectively, for γ = AR = 1 (circular vents), and, for fissure vents, γ decreases as AR increases.
For simplicity, the effects of topography were neglected, assuming the volcanic vent to be located on a flat surface. Furthermore, we assumed that the solid pyroclasts were dynamically and thermodynamically well coupled with the gas phases (i.e., volcanic gas and entrained air), ignoring particle sedimentation and decoupling (see Cerminara, Esposti Ongaro, & Berselli, 2016, Cerminara, Esposti Ongaro, & Neri, 2016. In the presented models, the transport equations that describe the mass, momentum, and energy conservation as well as the equations of state were solved numerically. The numerical pseudo gas model was based on the time-dependent solution of the compressible flow Navier-Stokes equations. The model used implicit large-eddy simulation (ILES) approach where the inherent numerical dissipation act as an implicit subgrid model. We confirmed that when spatial resolution is sufficiently high, both the numerical results with and without the LES approach correctly reproduce the spreading rate of jets observed in the experiments (Suzuki et al., 2005). The simulation results, obtained by the present model, were compared with the other models that used the LES approach (see also the intercomparison of 3D models in Suzuki, Costa, Cerminara, et al., 2016). For the detail of the model, see Suzuki et al. (2005) and Suzuki and Koyaguchi (2013).
The computational domain had a horizontal extent of~100 km and a vertical extent of~60 km. Non-uniform structured grids covered the physical domain. The minimum grid cell size was used for the vent and ground surface; it was smaller than L min /6, where L min is the minimum length scale of the vent (Table 1). The grid cell size increases with increasing horizontal distance from the crater and with increasing vertical distance from the ground surface. The grid cell size was increased at constant rates of 0.005 and 0.02 units horizontally and vertically, respectively.
A representative tropical atmosphere was used for the initial atmospheric conditions. Atmospheric density and pressure were stratified for a given temperature gradient, and, for simplicity, we assumed calm (i.e., no wind) conditions. We applied no-slip conditions at the ground boundary and free outflow/inflow conditions at the other boundaries. We focused on two eruption intensities with parameters close to the critical condition required for column collapse (Costa et al., 2018). The first eruption type had a mass flow rate (MFR) of 10 9 kg s −1 , and the second had an MFR of 10 9.5 kg s −1 . The vent exit velocity, temperature, water fraction, and density were set at 256 m s −1 , 1,053 K, 0.06, and 3.5 kg m −3 , respectively (see Costa et al., 2018). Circular and fissure (i.e., rectangular) vent geometries of different dimensions were considered, and the length of the fissure was varied (3.4, 10, and 15 km for MFR = 10 9.0 kg s −1 ; 10 and 20 km for MFR = 10 9.5 kg s −1 ). For fixed vent geometry, exit velocity, and magmatic properties, the vent radii and fissure widths can be determined for a given MFR. The input parameters are listed in Table 1.

Simulation Results
The simulation results show that the flow patterns of volcanic plumes are controlled by the volcanic vent geometry. Here, we describe the simulation results as a function of the AR (see Table 1).
For the magmatic properties considered in this study, when the MFR was set at 10 9 kg s −1 , the eruption clouds developed stable columns and did not produce clearly recognizable pyroclastic flows (Figure 1), regardless of the vent shape. The plumes from circular vents (γ = AR = 1) were associated with highly concentrated flow of the erupted mixture above the vent (at a height about four times the vent diameter), termed "radially suspended flow". When the vent radius is large, the erosion by the outer shear layer of the flow does not reach the central axis at the height where the initial momentum is exhausted (defined as H mmt ). The dense unmixed cloud stops to rise and spreads radially at this height, resulting in the formation of the RSF. The large vortices at the top of the RSF entrain a significant amount of ambient air, and the mass fraction of the erupted mixture was significantly reduced in this region (see supporting information, Movies 1a-1c). For a 3.4-km-long fissure vent (AR = 10), the eruption cloud also exhibited RSF structures (Figure 1b). For high-aspect-ratio fissure vents (AR = 91 and 200), a sheet-like buoyant flow (i.e., a quasi-2D vertical flow) developed, and no RSF was observed (Figures 1c and 1d; Movies 3a-3c and 4a-4c). In these latter cases, the mass fraction of the erupted mixture gradually decreased with height in contrast to the sudden decrease at the top of RSF. Among our simulations, the cases having AR from 1 (circular vent) to 22 belong to the RSF regime, while the cases with high aspect ratios (AR = 91 or 200) belong to the non-RSF regime.
Changes in plume height over time can reflect the type of flow pattern (Figure 2a). In Runs 1-4 (see Table 1), the height of the plume increased until~400 s. This period is considered to represent the initial stage of plume development. In the RSF regime (Runs 1 and 2), the speed at which the plume increased in height slowed at~5-10 km at~100 s but subsequently recovered to the initial rising speeds. These observations Note. γ = the deviation of the vent shape from a true circle; ΔX max = maximum grid size; ΔX min = minimum grid size; AR = the aspect ratio; L max = maximum length scale; L min = minimum length scale; MFR = mass flow rate.  Figure 2b illustrates the vertical distribution of mass fraction of the erupted mixture for the series of simulations with MFR = 10 9 kg s −1 . In the RSF regime (Runs 1 and 2), the average mass fraction of the erupted mixture was high near the vent (0.9-1.0 in Run 1, 0.8-0.9 in Run 2). Above the RSF, the mass fraction significantly decreased. The mass concentration in the buoyant region increased as AR increased (see also Figures 1a and 1b). In the non-RSF regime (Runs 3 and 4), the erosion by the outer shear layer reaches the central axis, and the unmixed core disappears at a height (defined as H core ) less than H mmt . In this regime the mass concentration decreased gradually with increasing height, but the average mass fraction within the buoyant region remained high (~0.6 for Run 3,~0.3 for Run 4). In contrast to the RSF regime, the average mass concentration in the buoyant region was observed to decrease as AR increases (Figure 2b). The observation that plume height increased with increasing mass concentration within the buoyant plume ( Figure 2 and Table 2) is consistent with previous results from 1D models (e.g., Morton et al., 1956). , whereas the dashed ones represent those for 10 9.5 kg s −1 . The simulation results of Run 1hr (high resolution) with high spatial resolutions show qualitatively similar features as that of Run 1. The time window used to compute averages in (b) is 600-800 s. The black curve represents the result of the 1D model (Woods, 1988). Note. Ri(L 0 ), Ri(L min ) = Richardson number based on L 0 and L min , respectively. H max is the maximum plume height estimated from time average for 600-800 s in Figure 2a. H core was estimated from the highest level where the mass fraction of the erupted mixture is 0.9 in Figure S3.
When the MFR is set at 10 9.5 kg s −1 , the volcanic plumes were less stable than those with an MFR value of 10 9 kg s −1 (Figure 1e). For the circular vent scenario (Run 5), a small portion of the ejected mixture was seen to collapse from the top of the RSF, falling to the ground, but the majority of the ejected mixture became buoyant (Figure 1e; Movies 5a-5c). This scenario is indicative of an incipient collapsing regime whereby the collapsing mixture spread on the ground, entrained ambient air, and rose upward where it coalesced with the main buoyant plume (see Costa et al., 2018). For the 10-km-long fissure (Run 6; AR = 28), the simulated cloud also showed a partial column collapse from the top of the RSF (Figure 1f; Movies 6a-6c). In contrast, the volcanic plumes from vents with higher aspect ratios (i.e., the 20-km-long fissure, Run 7; AR = 111) were more stable and did not develop RSF structures or pyroclastic flows (Figure 1g; Movies 7a-7c). In this scenario, the plume showed a gradual decrease in mass fraction with increasing height, whereas the plumes in the RSF regime showed the sudden decrease (Figure 2b).

Discussion
In this section, we discuss how vent geometry affects the dynamics of volcanic plumes. Our simulation results indicate that the aspect ratio of the vent (AR) significantly affects volcanic plume flow patterns, which in turn controls maximum height and plume stability.
One of the most important controls on the dynamics of volcanic plumes is how the erupted mixture entrains ambient air. The amount of entrained air depends on the entrainment velocity and the peripheral surface area of the plume; the latter directly depends on vent geometry only near the vent, say below the buoyant region. Morton et al. (1956) proposed that the entrainment velocity is proportional to the velocity of a turbulent jet or plume in a uniform environment (i.e., entrainment hypothesis). Some studies have modeled fissure eruptions (Glaze et al., 2011(Glaze et al., , 2017Stothers, 1989;Woods, 1993) based on the entrainment hypothesis assuming that the proportionality constant (i.e., entrainment coefficient) itself is independent of vent geometry. In this case, the amount of entrained air depends only on the entrainment area (hence on γ). A small value of γ represents a large entrainment area and γ = 1 represents the minimum value of the entrainment area. As mentioned above, γ decreases with increasing AR (see Table 1).
The fact that for the non-RSF regime, averaged mass concentration in the buoyant region decreases as AR increases (Figure 2b), leading to a decrease in the maximum plume height (Figure 2a) is consistent with the above assumption that the amount of entrained air below the buoyant region increases as the entrainment area increases. In the RSF regime, however, the mass concentration increases as AR increases (Figures 1a, 1b, and 2b) even though the entrainment area increases (i.e., γ decreases). This is the opposite effect of that observed in the non-RSF regime, suggesting that the entrainment coefficient itself is strongly dependent on the flow pattern and, particularly, the presence of RSF (cf. . The critical condition that separates the RSF and non-RSF regimes is determined by the flow geometry near the vent. This condition is the balance between two length scales: (1) the height where the unmixed core disappears, H core , and (2) the level where the initial momentum is exhausted, H mmt Suzuki & Koyaguchi, 2012). When H core < H mmt , the highly concentrated core disappears before the initial momentum is exhausted, and the flow is associated with a non-RSF regime. On the other hand, when H core > H mmt , the highly concentrated core develops RSF at H mmt . It is important to highlight that H mmt is a function of velocity and density of the erupted mixture (i.e., magmatic properties) at the vent and is independent of vent geometry. For our simulations, H mmt is estimated to be 5 km ( Figure S4). In contrast, H core is strongly dependent on vent geometry and is proportional to the vent radius for a circular vent or the minimum length scale of the vent (L min ) for a fissure; the value of H core /H mmt decreases as AR increases (see Figure S3 and the value of H core in Table 2). As a result, for the high-aspect-ratio scenario (Run 7 with AR = 111), H core /H mmt is so small that a well-mixed buoyant plume, belonging to the non-RSF regime, is generated.
Our results also suggest that column collapse depends on whether the plume belongs to the RSF or non-RSF regime. According to previous studies (e.g., Woods, 1988), column collapse is controlled by the balance of momentum and buoyancy fluxes at the source (i.e., the source Richardson number), and, when the magma properties and exit velocity are fixed as in the present cases, its critical condition is primarily determined by the amount of air entrained in the plume. In the RSF regime, in which a highly concentrated core remains up to H mmt (Runs 1 and 2; see Figure 2b), column collapse is determined by the amount of air entrained within the RSF. In this regime, as the amount of entrained air in the RSF decreases as AR increases, the plume is expected to be less stable with increasing AR. As AR further increases, a transition from an RSF to a non-RSF regime occurs. In the non-RSF regime (Runs 5-7; Figures 1e-1g), the erupted mixture around the central axis of the vent is well mixed with ambient air, resulting in a steady decrease in density with increasing height, which enhances the stability of the volcanic plume. Unlike the case of the RSF regime, the plume is expected to be more stable with increasing AR within this regime. The observed transition from an RSF to non-RSF regime with increasing AR provides evidence for a newly identified control of vent geometry on the requirements for column collapse.
The above features associated with the RSF structures cannot be described by 1D models based on the entrainment hypothesis (Morton et al., 1956), since the flow pattern around the RSF is highly complex with a multidimensional structure (Cerminara, Esposti Ongaro, & Neri, 2016;. However, 1D models are still useful for predicting the stability and height of volcanic plumes if the critical condition between the RSF and non-RSF regimes and the changes in the entrainment efficiency associated with the transition between the two regimes are considered. Suzuki and Koyaguchi (2012) showed that, for circular vents, the critical condition separating the RSF and non-RSF regimes is given by a formula that uses the source Richardson number (equation 6 in their paper). To apply this formula to fissure eruptions, the effect caused by reduced H core /H mmt due to the minimum length of the vent, L min , should be taken into account; in practice, the vent radius used in the calculation of Richardson number should be replaced with L min (Table 2 ). Once this Richardson number is introduced, we can discuss the critical condition between the RSF and non-RSF regimes for a wide range of MFR. Our study reveals the opposite dependency of plume height on vent aspect ratio between the RSF and non-RSF regimes. To calculate plume height, future studies need to estimate the effective value of the entrainment coefficient for both the RSF and non-RSF regimes.
The critical condition separating the RSF and non-RSF regimes and that for column collapse can be affected by other factors. The particle decoupling from gas phase can alter the entrainment properties of plumes (e.g., Cerminara, Esposti Ongaro, & Neri, 2016;Jessop & Jellinek, 2014;Jessop et al., 2016). The crater geometry, such as the opening angle of crater, controls the decompression-compression processes and the exit velocity at the crater top (Carcano et al., 2014;Cigala et al., 2017;Kieffer, 1989;Ogden et al., 2008;Woods & Bower, 1995) and the entrainment properties (e.g., Cerminara, Esposti Ongaro, & Neri, 2016;Jessop & Jellinek, 2014). Although these factors may quantitatively change the critical conditions of RSF/non-RSF and stable/collapse regimes, the present results that column collapse depends on whether the plume belongs to the RFS or non-RSF regime as function of the vent shape is considered to be robust because the difference of these regimes is caused by the fundamental features that the clouds are ejected from finite length-scale vent and their buoyancy changes from negative to positive. Finally, further studies should consider how plume dynamics change when magma properties (e.g., temperature, water content) and vent geometries are varied together.

Conclusions
Our study showed that vent geometry has two different effects on entrainment process: the peripheral surface of the plume and the flow structure (e.g., the formation of RSF). The combination of these two effects have a fundamental control on the stability and height of volcanic plumes. Volcanic plumes of high mass flow rate from circular or low-aspect-ratio vents are associated with RSF in the lower parts of the plume, near the vent, and are less stable than plumes from high-aspect-ratio vents. For the RSF-type plume, the amount of air within the plume is controlled mainly by the entrainment efficiency around the RSF. The entrainment efficiency decreases and the plume height increases as the aspect ratio increases. In contrast, the amount of ambient air entrained into non-RSF plumes depends on the entrainment area. In this regime, the entrainment area increases and the plume height decreases as the aspect ratio increases. Column collapse is also critically dependent on whether or not the plumes are associated with RSF.