Multiphase turbulent flow explains lightning rings in volcanic plumes

The January 15, 2022 HTHH eruption produced an extraordinary 57–58 km-tall plume, an umbrella cloud expanding to ~400 km diameter in under an hour, and globally propagating geophysical waves. It also generated unprecedented lightning activity (~400,000 strokes over 6 hours) that formed concentric ring patterns centered on the vent. These regular, periodic lightning rings differ from gaps seen in supercell thunderstorms and raise the question of their origin. The study investigates whether turbulence in a buoyant, particle-laden plume within a stably stratified atmosphere can produce the observed lightning rings via particle clustering, and whether lightning ring properties can serve as proxies for otherwise obscured eruption dynamics and parameters.

Background on HTHH activity notes prior Surtseyan eruptions and escalating activity leading to the main 15 January event. Volcanic electrification can arise from size-dependent triboelectrification of colliding ash particles, collisions of hydrometeors formed in situ, and processes associated with magma-water interaction during submarine eruptions. Lightning rates can reflect particle acceleration, turbulence, and plume extent, but the role of turbulence in electrification needs clarification. Turbulent flows can preferentially concentrate inertial particles, enhancing collisions relative to homogeneous distributions. Classical plume theory (Morton, Taylor & Turner, 1956) treats buoyancy-driven plumes/thermals with entrainment and Boussinesq assumptions, predicting scaling for plume radius and maximum height, and umbrella cloud spread as gravity currents. Volcanic eruptions add complexities (gas-thrust jet region, winds, external water entrainment). External water can increase or decrease buoyancy depending on entrainment fraction. Prior 3D multiphase simulations included gas-particle decoupling and compressibility, and noted turbulence-driven clustering but focused on large-scale plume dynamics. Minimal models have been advocated to isolate essential mechanisms without subgrid parameterizations.

The study uses a minimal 3D numerical model of a buoyant plume in a dry, linearly stratified, incompressible Boussinesq fluid with passive, heavy inertial point particles under linear Stokes drag. Flow equations (Boussinesq momentum and density fluctuation with linear stratification) are solved in a periodic domain L×L×H = 188×188×94 km³ using a high-order pseudo-spectral code (GHOST) at 1024×1024×512 resolution (with convergence checks at lower resolutions). A Gaussian ellipsoidal source of buoyancy (light/hot fluid) near the bottom provides a constant energy input, generating an updraft plume that overshoots the level of neutral buoyancy (LNB) and spreads to form an umbrella cloud. Background Brunt-Väisälä frequency is N = 0.01 s⁻¹ (gravity-wave period ~100 s). Updraft velocities at the plume center reach ~150 m s⁻¹; the domain-wide rms turbulent velocity at late times is U ≈ 26 m s⁻¹, while in the umbrella cloud U_rms ≈ 61 m s⁻¹. The Froude number Fr = U/(LN) ≈ 0.01 (for L = 188 km). Kolmogorov scales (set by chosen ν, κ with turbulent Prandtl ~1): τη ≈ 35 s, η ≈ 115 m, uη ≈ 3.2 m s⁻¹; ν ≈ 378 m² s⁻¹. Energy input rate per unit mass ε ≈ 1.2 W kg⁻¹ gives total É ≈ 1.6×10¹² W (Gaussian source volume characterized by vertical dispersion 2.8 km and horizontal 5.6 km), comparable to a mass discharge ~2×10⁸ kg s⁻¹ with 800 K excess temperature and heat capacity 1000 J kg⁻¹ K⁻¹. Particles obey the heavy point-particle limit of the Maxey–Riley equation with linear drag; gravity and settling are neglected in baseline runs and included for one case. Four particle sets (3 million particles each) are integrated: species A (St=1), B (St=0.1), C (St=0.01) initially in a thin basal layer at rest, and initially suspended particles (ISP, St≈1 like B parameters) homogeneously distributed between z=10–60 km. Stokes number St = τp/τη; typical τp: 35 s (A), 3.5 s (B), 0.35 s (C). Particle sizes are set to keep realistic ratios to Kolmogorov scales (d_p/η ≲ O(1)) and Re_p<1, ensuring validity of linear drag. Diagnostics include particle number densities n_i(r) vs radius r, squared vorticity ω²(r) as a turbulence proxy, and collision proxies N_ij(r) ∝ n_i n_j (v_i v_j)^1/2 between species pairs (e.g., A–C, ISP–C). Temporal evolution of ring positions and cloud edge are tracked; gravity-wave and turbulence influences are assessed via particle velocity spectra in umbrella regions. Visualizations use VAPOR.

• Turbulence-induced inertial clustering produces strong particle accumulation in two regions: the central column (r ≲ 20 km) and an annular ring at r ≈ 40 km across particle sizes (St = 1, 0.1, 0.01) and for initially suspended particles. This creates an annular gap between core and ring, matching observed lightning ring morphology.
• Radial profiles at t = 43 min show all bottom-injected species (A, B, C) have density maxima at r = 0 and at r ≈ 40 km, coincident with peaks in ω²(r). ISP also peak at r ≈ 25 and 40 km and tend to homogeneous values at large r.
• Collision proxies N_AC and N_ISP,C peak near the axis and at r ≈ 40 km, indicating enhanced probability of inter-size collisions at the ring, favorable for electrification.
• Ring behavior over time: while the umbrella cloud edge continues to expand, the ring repeatedly forms and fluctuates around a mean radius ≈ 40 km, with linear outward propagation episodes (speed ≈ 65 m s⁻¹) similar to observations. Multiple ring occurrences are recreated spontaneously even with constant source flux.
• Expansion scalings: initial approximately linear expansion (∝ t) transitions to slower growth consistent with gravity-current theory for umbrella clouds (∝ t^0.7, later ∝ t^0.4). These trends are consistent with observations from HTHH.
• Gravity-wave vs turbulence effects: Particle velocity spectra exhibit a peak near N (Brunt-Väisälä frequency) in vertical motions, indicating vertical confinement/oscillation by gravity waves (~100 s periods), whereas horizontal motions show Kolmogorov-like spectra with lower characteristic frequencies, implying turbulence dominates horizontal ring displacement.
• Gravity/settling: Including gravity for St=1 still yields a density ring at r ≈ 40 km and can increase clustering in strongly turbulent regions; for St ≤ 0.1, settling speeds are smaller than turbulent velocities, so gravity has minor impact on clustering patterns.
• Observational linkage: The model reproduces key features of HTHH lightning rings (ring radius ~40 km, expansion speed ~65 m s⁻¹, repeated ring formation). Analysis of spatiotemporal lightning density suggests significant explosive pulses around 04:14, 04:51, 05:34, and 08:33 UTC, consistent with independent seismic, tsunami, and acoustic observations.
• Quantitative context: HTHH lightning peaked at ~400,000 strokes over 6 h (>5000 strokes/min), with rings concentric about the vent; the simulations relate ring locations to turbulence maxima and particle clustering zones.

The simulations demonstrate that multiphase turbulence in a buoyant, stratified volcanic plume can preferentially concentrate particles in a central updraft core and in an annulus where turbulence intensity peaks, naturally producing lightning ring structures and an intervening lightning gap. This mechanism differs from thunderstorm 'lightning holes,' which track drifting updraft surges; in HTHH, the central region over the vent remained lightning-active while rings appeared and dissolved in a persistent concentric band. Enhanced collisions between differently sized particles (e.g., ash and hydrometeors/ice) are expected where densities co-locate, supporting non-inductive electrification. Stratification is essential: it confines motion vertically and helps establish the annular high-vorticity region at a relatively fixed radius, while turbulence drives the horizontal organization and repeated ring formation. The modeled expansion of both the umbrella cloud and the lightning ring, and the ring’s fluctuation around ~40 km, provide a physical basis to interpret lightning observations as proxies for in-plume dynamics. Importantly, ring expansions can arise from intrinsic plume-turbulence fluctuations even under constant source flux; thus, distinguishing rings originating at the vent (indicative of new explosions) from those initiating at finite radius (intrinsic variability) allows inference of eruption pulses obscured by the umbrella.

A minimal, incompressible Boussinesq plume model with heavy inertial particles reproduces the formation, evolution, and persistence of lightning rings observed during the HTHH eruption. Turbulence-driven particle clustering in a stratified environment concentrates particles in a central core and an annular region (~40 km radius), maximizing inter-size collision probabilities and enabling strong electrification there. The model captures umbrella cloud radial growth and ring oscillations, showing that repeated ring expansion/contraction can occur without unsteady source flux and allowing identification of eruption pulses despite optical obscuration. The work proposes using lightning ring location, size, speed, and persistence as proxies for plume turbulence and eruption parameters. Future efforts should incorporate more realistic physics—moist processes and phase changes, atmospheric profiles, compressibility and supersonic vent jets, two-way particle-fluid coupling, particle interactions and charging, and time-varying source conditions—to enable quantitative retrieval of eruption source parameters from joint lightning and cloud observations.

The model is intentionally simplified: incompressible Boussinesq framework with linear stratification; dry plume without condensation, hydrometeor formation, or phase changes; one-way coupling with non-interacting point particles (no particle-particle collisions or feedback on the flow); predominantly neglects gravity/settling (included only in one case); constant buoyancy source after onset (no explicit pulsatory vent forcing); no gas compressibility or supersonic jet dynamics in the gas-thrust region; and significantly inflated Kolmogorov scales (η, τη) due to computational constraints compared to real eruptions. These limitations affect direct quantitative comparisons to observations and preclude explicit electrification/charging calculations.