Numerical simulations of the latest caldera-forming eruption of Okmok volcano, Alaska

The 2050 ± 50 ¹⁴C yBP caldera-forming eruption of Okmok volcano, Alaska, had a global atmospheric impact with tephra deposits found in distant Arctic ice cores and a sulfate signal found in both Greenland and Antarctic ice cores. The associated global climate cooling was driven by the amount of sulfur injected into the stratosphere during the climactic phase of the eruption. This phase was dominated by pyroclastic density currents, which have complex emplacement dynamics precluding direct estimates of the sulfur stratospheric load. We simulated the dynamics of the climactic phase with the two-phase flow model MFIX-TFM under axisymmetric conditions with several combinations of mass eruption rate, jet water content, vent size, particle size and density, topography, and emission duration. Results suggest that a steady mass eruption rate of 1.2–3.9 × 10¹¹ kg/s is consistent with field observations. Minimal stratospheric injections occur in pulses issued from the central plume initially rising above the caldera center and from successive phoenix ash-clouds caused by the encounter of the pyroclastic density currents with topography. Most of the volcanic gas is injected into the stratosphere by the buoyant liftoff of dilute parts of the currents at the end of the eruption. Overall, 58–64 wt% of the total amount of gas emitted reaches the stratosphere. A fluctuating emission rate or an efficient final liftoff due to seawater interaction is unlikely to have increased this loading. Combined with petrological estimates of the degassed S, our results suggest that the eruption injected 11–20 Tg S into the stratosphere, consistent with the subsequent climate response and Greenland ice sheet deposition. Our results also show that the combination of the source Richardson number and the mass eruption rate is able to characterize the buoyant–collapse transition at Okmok. We extended this result to 141 runs from 10 published numerical studies of eruptive jets and found that this regime diagram is able to capture the first-order layout of the buoyant–collapse transition in all studies except one. An existing multivariate criterion yields the best predictions of this regime transition.