Submarine terraced deposits linked to periodic collapse of caldera-forming eruption columns

Catastrophic, caldera-forming explosive eruptions generate hazardous ash fall, pyroclastic density currents and, in some cases, tsunamis, yet their dynamics are still poorly understood. Here we use scaled analogue experiments and spectral analysis of well-preserved concentric terracing of seafloor deposits built by submarine caldera-forming explosive eruptions to provide insights into the dynamics governing these eruptions and the resultant hazards. We show that powerful submarine eruption columns in collapsing regimes deliver material to the sea surface and seabed in periodic annular sedimentation waves. Depending on the period between successive waves, which becomes shorter with decreasing jet strength, their impact and spread at the sea surface and/or seabed can excite tsunamis, drive radial pyroclastic density currents and build concentric terraces with a wavelength that decreases with distance, or deposits that thin monotonically. Whereas the Sumisu (Izu–Bonin arc) caldera deposit architecture is explained by either a subaerial or deep-water model involving no interaction between sedimentation waves and the sea surface, those of the Macauley (Kermadec arc) and Santorini (Hellenic arc) calderas are consistent with a shallow-water model with extensive sedimentation wave–sea surface–seabed interactions. Our findings enable an explicit classification of submarine caldera-forming explosive eruption dynamics and quantitative estimates of eruption rates from their terraced deposits. Submarine terraced deposits of some caldera-forming explosive eruptions result from periodic collapses of the eruption column and can be used to estimate their source eruption rate, according to an analysis of such terraces and analogue experiments.

. In this article, we focus on the dynamics governing the delivery of dense mixtures of coarse pyroclastic material in SWs to the water surface (free surface) and seafloor, which forms the major volume of the erupted material 15,[20][21][22] . Our experiments do not model the delivery of fine ash to the upper troposphere and stratosphere by buoyant plumes of interstitial gas that escape the collapsing part of the eruption column where it stalls at or near the maximum 'fountain height', ~H ftn (refs. 26,27). We also neglect complex effects of transient, spatially varying layers of floating pumice (pumice rafts) on SW-free-surface interactions 28 (further discussion in Methods).
In 'deep-water' experiments, the water-layer depth h w ≫ H ftn (equation (5); Figs. 1e and 3a), whereas in shallow-water experiments h w ≤ H ftn (Figs. 1f and 3b and Supplementary Videos 1 and 2). Our previous work shows that in collapsing regimes in air or deep water, eruption columns collapse periodically as annular SWs to deliver material to an impact zone (Figs. 1e,g and 3a; fig. 20 in ref. 19). The excitation of SWs is driven by the periodic accumulation and release of dense material at H ftn , which excites oscillations of the fountain top at distinct fountain and environmental stratification frequencies 29 (equations (6) and (7), respectively). Usefully, the fountain frequency predicts the timing between SWs and resulting PDCs, which govern the architecture of the terraced deposits (Methods and Fig. 3a; fig. 17 in ref. 19).
For shallow-water experiment #11, four distinct oscillation periods are identified in the subaerial fountain-top height time series (Fig. 3c,d and Supplementary Fig. 2) and agree with predicted 1/f ftn , 1/N wtr and 1/N air from source and environmental conditions, as well as 1/f sw e measured by tracking SW descent speed and radius. Overall, we find that shallow-water fountains input energy to the stratification at f wtr ftn and excite the higher frequencies N wtr , N air , f air ftn and f sw e . Energy is dissipated at frequencies above f sw e predominantly through molecular mixing across the air-water density interface 30 .
In the absence of a water layer or under deep-water conditions, SWs carry a momentum flux that increases in proportion to H ftn . Depending on H ftn and the extent to which entrainment and mixing reduce the density excess driving the descent of these mixtures, SWs can either impact the ground or seabed as jets that drive intense scouring and erosion or become largely dissipated descending clouds of dense particles. Deposition from such relatively quiescent particle clouds is mostly by individual particle settling 19,31 . By contrast, where strong impacts occur, SWs evolve into axisymmetric PDCs that deposit regularly spaced and axisymmetric terraces (Figs. 1c,e,g and 2e,f). Geometrically, terraces are confined within a characteristic impact zone with radius r im from the jet margin that is comparable to the jet diameter.
Under shallow-water conditions, SWs descending through air impact the free surface and excite water waves (analogue tsunamis) (Fig. 1f,h and Supplementary Videos 1 and 2). The momentum-driven spread of the mixture and its overshoot across the free surface as an erosive jet can scour the seabed, depending on the water-layer height h w . A useful metric to characterize SW interactions with the seabed is, thus, the dimensionless penetration depth where H = H ftn for subaerial or deep-water conditions or H = h w for shallow-water conditions and l sw is a characteristic scale for the overshoot depth of an SW (equation (9)), which can also be inferred observationally. If D sw > 1, pyroclastic particle settling within the impact zone will form a nearly planar terrace. By contrast, if D sw ≤ 1, descending SWs impact and spread at the seabed within the impact zone as erosive jets that deposit backward-facing terraces. In the special limit D sw ≪ 1, intensive scouring produces a distinctive concave backward-facing terrace architecture (Fig. 4). In this case, D sw predicts also a 'scouring gravity currents such as landslides and PDCs be expressed so similarly at calderas with disparate structures, mechanical properties, slopes and slope stabilities and with distinct volcanic source geometries and eruptive histories? Recent analogue experiments on subaerial eruptions suggest that concentric terracing (Fig. 1c,d) is an intrinsic feature of the particle-fluid (multiphase) interactions governing the rise and descent of material in eruption columns, and spread along the ground thereafter, evolving in 'partial collapse' or 'total collapse' regimes ('collapsing regimes' hereafter) 19 . Most well-documented large CCF events in the geological record, including CCF eruptions in shallow-water environments, involve eruptive phases where half or more of the total erupted mass was delivered to PDCs 15,[20][21][22] . In collapsing regimes that produce terraced deposits, as much as 90% of material ejected into the atmosphere is delivered back to the ground through the excitation and descent of periodic, annular sedimentation waves 19 (SWs; Fig. 1e). Critically, the minimum radius from the vent of terraces built by this mechanism is predicted to be similar to the erupting jet radius taken at its mean rise height, H ftn . While this expectation is potentially consistent with the terracing observed ~1 km adjacent to the caldera rim (in the 'near field') at Sumisu caldera ( Fig. 1a), it is inconsistent with the relatively much larger radii of near-field terraces observed at Santorini (Fig. 1b) and Macauley calderas 13 ( Supplementary Fig. 1).
A key difference to subaerial eruption columns in collapsing regimes, however, is that the dynamics of submarine CCF events are modified by interactions with water layers, which cause the mixture to spread as it penetrates and descends 23,24 . In this article, we use analogue experiments on submarine eruptions through water layers of varying depth to show that terrace formation similarly occurs through the excitation of periodic SWs (Fig. 1c-h and Methods) and that this process is probably a generic feature of CCF events. We show that this process and its expression in the architectures of resulting deposits are highly sensitive to, and diagnostic of, both eruptive source and water depth conditions. We also show that the structural characteristics of terraces closest to the caldera rim record the effects of erosion and sedimentation related to the impact and spread of SWs at the seabed as PDCs, consistent with earlier studies 13 . Taken together, our results strengthen an emerging framework for classifying eruption column collapse deposits, quantitatively constraining eruption source parameters from qualitative observations of deposit architectures and understanding the timing and intensity of tsunami and PDC hazards related to CCF eruptions 19,25 .

Periodic sedimentation waves build terraces
CCF events at the Santorini, Sumisu, and Macauley (Supplementary Information) calderas 13 each delivered mixtures of ash, gases and water vapour to the atmosphere through a shallow water layer. Each also produced deposits marked by sharply defined, approximately concentric and periodic terraces with low-slope, backward-facing profiles (Figs. 1a,b and 2a,b and Supplementary Fig. 1). Terraces at Santorini and Macauley decline monotonically in wavelength from order 1,000 to 100 m over 10-16 km radial profiles. This pattern of wavelength decline is also evident in profiles from the northeastern terraces around Sumisu 11 (Figs. 1a and 2b). However, significant differences in the deposit architectures occur. Whereas terraces are approximately planar ~1 km from the caldera rim in the near field at Sumisu, near-field terraces with concave backward profiles at Santorini and Macauley calderas occur ~3−6 km from the caldera rim (Figs. 1a,b and 2a,b and Supplementary  Fig. 1). Furthermore, there is a relatively wide and deep trough between the rim of Macauley caldera and the start of the near-field terraces that is not present at Santorini or Sumisu ( Supplementary Fig. 1).
To investigate links among eruption column collapse dynamics, SWs and the construction of terraced deposits during CCF events in submarine settings, we conduct carefully scaled analogue experiments on turbulent particle-water fountains injected into water layers of varying depths and analyse the fountain-top height oscillations and Article https://doi.org/10.1038/s41561-023-01160-z radius' R sc = r sc /r 0 , where r sc is the radius of the scoured region of the deposit and r 0 is the vent radius (Fig. 1c). Of practical value (SWs and shallow water layers in collision), field measurements of the scouring radius can constrain the water depth and momentum flux carried by SWs. Furthermore, if the frequency at which SWs impact the free surface is higher than the frequency at which spreading SWs descend to     11 . b, Bathymetry of the shallow submarine terraced pyroclastic deposit surrounding Santorini caldera (Greece, 3.6 thousand years ago 10,14 ). Positive and negative latitude values correspond to North and South of the equator, respectively. Similarly, positive and negative longitude values correspond to East and West of the prime meridian, respectively. c, Left image shows terraced deposit built by deep-water experiment #10 (Supplementary Table 1) with marked scouring radius r sc and terrace wavelengths (λ). Right figure shows a vertically exaggerated digital elevation model of deposits, where radial particle ridges are marked by white line. d, Shallow-water experiment #6 deposit marked similarly as in c. e,f, Image pairs of a descending SW in deep-water experiment #32 (e) and shallow-water experiment #6 (f) where r sw is the SW radius, l sw is the SW jet entrance length into the water layer (equation (9)) and t is time. g,h, Conceptual model cartoons of a descending SW in a deep-water (g) and a shallow-water (h) experiment. Panel a reproduced with permission from ref. 11, Springer Nature Limited.

SWs and shallow water layers in collision
A number of key geometric properties of terraces observed around Sumisu, Santorini and Macauley calderas place restrictive constraints on the mechanics of their formation. In plan view, broad, regular concentric terraces with scalloped edges have a wavelength that declines with distance from caldera rims ( Fig. 1). Furthermore, the occurrence of intermittent tsunamis is often associated with historic shallow-water caldera eruptions [32][33][34] , suggesting an intermittent tsunamigenic process during eruptions.
Our analogue terraced deposits are diagnostic of eruption intensity in collapsing regimes and have planforms and radial profiles that are sensitive to the water-layer and SW penetration depths D sw . Consistent generally with deposits at Sumisu, Santorini and Macauley calderas, analogue deposits are constructed of regularly spaced, concentric and backward-facing terraces with broad, low-slope top surfaces, sharp crests and relatively short, high-slope downstream sides. Evident visually (Figs. 1c,d, 4 and 5) and spectrally (Fig. 2e,f), periodic terrace widths also decline with radial distance from the near to far field of experimental deposits. Akin to deep-water experiments under intermediate or strong eruptive conditions, the backward-facing Sumisu deposit terraces are also characterized by azimuthal scalloping along terrace edges over scales smaller than terrace widths in the near field and comparable to the terrace widths in the far field. By contrast, concavity in near-field Santorini and particularly Macauley terraces signal both h w ≪ H ftn and D sw < 1 conditions and require that these events be very strong eruptions (Figs. 1b, 2a and 5). Indeed, the ~2-km-wide concavity containing eruption-derived material at Macauley is consistent with intense excavation within an impact zone under potentially extreme D sw ≪ 1 conditions. Quantitatively, water depth and eruption fountain strength are indicated from the deposit architecture through a Richardson numberparticle volume fraction (-Ri 0 ↔ ϕ 0 ) regime diagram derived from our experiments ( Fig. 5 and equations (2) and (3) in Methods). In deep-water regimes, the number and maximum radii of terraces decrease as fountains become weaker or more concentrated. In shallow-water cases, the water depth exerts the primary control over the extent to which SWs scour or deposit sediment in the near-field deposit (Fig. 4). Lateral scouring and concavity are intensified as D sw → 0, whereas near-field terraces become increasingly planar as D sw approaches and exceeds 1. Deposit architecture can be combined with plausible estimates of eruption source parameters to constrain the fountain strength, particle volume fraction and mass eruption rate (MER) of CCF eruptions ( Fig. 5 and Methods). MER exerts a profound control over eruptive behaviour and is the greatest source of uncertainty in computational models of   Spectral results of deposit profiles shown in Fig. 1a,b show a remarkable radial symmetry with power concentrated in the larger near-field terraces and a fall off with distance that is relatively smooth at Sumisu. e,f, Power spectra of deposit topography profiles shown in Fig. 1c,d for shallow-water experiment (e) and deep-water experiment (f) with near-and far-field terrace wavelength peaks marked.
Article https://doi.org/10.1038/s41561-023-01160-z eruptions aimed at understanding PDC hazards as well as volcanic effects on climate 24,35,36 . To our knowledge this is the first attempt at a reconstruction of CCF eruption source parameters from deposit architectures.

Generic model for terrace formation by SWs
We combine experimental results related to the periodic delivery of material in SWs with established models of flow and sedimentation by PDCs and turbidity currents to construct a generic conceptual model for the formation and evolution of terraced deposits related to successive SWs 19,37-39 (Fig. 6). During a shallow-water eruption, an initial SW impacts the free surface and spreads, generating interfacial waves (analogue tsunamis) (Fig. 6a). If D sw ≪ 1, the SW impacts and scours the seabed over a distance r sc (Fig. 4), deposits material to form a terrace ( Fig. 1) and drives a ground-hugging gravity current (Fig. 6a-d). Subsequent SWs impact and further scour the first terrace before they spread and deposit material downstream (Fig. 6e). Successive flows are modulated by interactions with evolving deposit topography through the production and draining of standing waves formed at hydraulic jumps at the downstream edges of terraces (Fig. 6e). A similar overall evolution occurs in deep-water experiments where D sw ≤ 1, albeit with less scouring of the deposit by successive SWs (deep-water experiments in Figs. 5 and 6g). It is important to note that this model applies only to deposition on relatively shallow and smooth slopes that are below the angle of repose of the particle size distribution carried by PDCs. The primary effect of a shallow water layer on this terrace formation process is to cause SWs to stall and spread as they overshoot the free surface before their descent and impact with the seabed. This process introduces an important additional link between deposit architecture and the eruption dynamics through a 'collapse frequency' number ℱ C (equation (10)). Key to the evolution summarized in Fig. 6 is the time between successive SWs. If ℱ C ≫ 1, SWs collapse in quick succession to combine or otherwise interact strongly before descending to the seabed or at the seabed, periodic terracing is inhibited and eruptions produce massive deposits, which is observed for the weakest events in our experimental results (Fig. 5). By contrast, as ℱ C → 0, the period between column collapses is long enough that PDCs can flow out of the SW impact zone r im as discrete events to build terraces in the near-field deposit with distinct bedding. In this case, however, if D sw → 0, deposit terraces adjacent to the caldera rim will be scoured (Figs. 4 and 5). Taken together, ℱ C and D sw provide explicit links among eruption source parameters, column dynamics, SWs, PDCs and deposit architecture.
Quantitatively, our results suggest that strong eruption columns in collapsing regimes will drive periodic rather than continuous PDCs to make terraced deposits where -Ri 0 ≳ 10 −4 and ℱ C < 1 (Fig. 5). Such dynamics involving intense SW-water surface interactions will lead also to tsunamis 32 and potentially to a greater likelihood for over-water and amphibious PDCs 40 . More generally, the inherent control of ℱ C over deposit architecture is not restricted to submarine events. Akin to the deep-water experiments, low-strength CCF eruption columns on land will also have a greater proclivity to produce massive deposits, consistent with previous studies 25, 41,42 . Whether subaerial terraced deposits     Fig. 3). Building on our previous study 19 , here we constrain minimum -Ri 0 − ϕ 0 values where terracing is expected to occur for subaerial and submarine eruption columns in collapsing regimes. Future studies should investigate the maximum -Ri 0 − ϕ 0 eruption column parameter space where ℱ C transitions from supplying discrete to continuous mass fluxes to spreading PDCs and whether this is expressed as a transition in near-field deposits from bedded and terraced architectures to massive architectures absent of terraces.
Our experimental results show that qualitative observations of terraced deposit architectures related to CCF events are linked to underlying eruption column dynamics and can quantitatively constrain source parameters of eruption columns in collapsing regimes. The potential to use these data to bound MERs, which are critical inputs for models and historically challenging to infer from observations 36 , provides exciting ways to study both submarine and subaerial eruption deposits. More broadly, understanding how the dynamics of PDCs driven by periodic SWs are potentially modified by evolving particle-fluid coupling regimes will also have implications for the flow of turbidity currents along the seafloor and rock and snow avalanches on land 43 .

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01160-z. a-e, Conceptual model cartoons of terraced and radial particle ridge deposit formation. Ground-hugging gravity currents shown in blue graded to light blue colours with particles shown in brown. The first three SWs are marked as SW 1 , SW 2 and SW 3 . f, Terrace modification by SWs and ground-hugging gravity currents at later times. g, SWs spreading as groundhugging gravity currents developing groovy instabilities 39 during a deep-water experiment with only fine particles (further discussion in Supplementary  Information). A surging ground-hugging gravity current is marked 'S', with lobe 'L' and cleft 'C' structures of other surging currents marked. Cross-section cartoon of streamwise helical flow corresponds to black boxes in d and g. Top view of an annular nozzle deposit from our previous work 22 showing scoured regions between terraces of decreasing wavelength with radial distance. Article https://doi.org/10.1038/s41561-023-01160-z

Methods
Here we describe our experimental method and scaling parameters that characterize eruptive and experimental source parameters 44 . We then show how the periodic properties of deposit terrace profiles are analysed to identify near-and far-field terrace wavelengths.
To model submarine explosive eruptions, we inject turbulent mixtures of water and silica particles with a constant flow rate into water layers deeper or shallower than the jet-rise height (deep-water and shallow-water experiments, respectively; Figs. 1e,f and 3a,b). For deep-water experiments, flows are injected into a density-stratified saltwater layer whereas for shallow-water experiments, they are injected into a constant-density water layer. We use 'fine' d p = 75 ± 25 μm and 'coarse' d p = 225 ± 25 μm diameter particles with densities of ρ f = 2,525 kg m -3 and ρ c = 2,693 kg m -3 , respectively, that simulate the inertial effects of fine and coarse ash on volcanic jet dynamics 22,45 . Applying experimental methodologies detailed in ref. 19, we characterize analogue eruption column behaviour quantitatively with a combination of high-resolution colour and greyscale images captured at 30 fps and 124 fps, respectively. Digital elevation models of analogue terraced deposits are produced after carefully draining experimental tanks with structure-from-motion photogrammetry carried out with Agisoft Metashape software.
The rise of both analogue and natural jets requires entrainment (and turbulent mixing). For entrainment to occur, work is extracted from the velocity field imparted at the mixture source to penetrate, overturn and engulf less-dense ambient fluid or atmosphere. The underlying balance between stabilizing buoyancy and driving inertial forces is captured by a source Richardson number 19,46,47 , Here, g ′ 0 = g(ρ a − ρ 0 )/ρ a is the buoyancy of the jet mixture where ρ a is the density of ambient water or atmosphere and ρ 0 is the bulk density of the mixture, r 0 is the jet radius and u 0 is the jet vertical speed, all taken at the source of the jet where it enters the ambient fluid. Usefully, −Ri 0 can be defined for multiphase jets at their source in a water layer, in an air layer or at the water-air interface where submarine jets rise through a shallow water layer to enter an air layer. In general, where -Ri 0 > 10 −4 , jets rise in collapsing regimes to a maximum height H ftn from which jet mixtures collapse periodically through the excitation of SWs 19,29,48 .
Our deep-water and shallow-water analogues address the key control of eruptive source conditions expressed through -Ri 0 . In deep-water/subaerial regimes, without considering additional thermal buoyancy fluxes arising through the heating of entrained atmosphere 49 and water-vapour condensation 50 , our predicted H ftn in terms of the source conditions alone is a lower bound. However, the magnitude of our underestimate is unclear. The extent to which condensation enters quantitatively is equivocal 51 . Furthermore, whether these thermal effects supply a sufficiently large enough buoyancy flux to the mixture to cause an order of magnitude or more change in the mass flux partitioned between spreading ash clouds and PDCs, or drive the rise of gas and fine ash out of the mixture, at or before jets rise to H ftn is also unclear.
For eruptions through shallow water layers with D SW ≤ 1, our experiments are strictly analogous to sonic or subsonic (pressure-balanced) eruptions where turbulent instabilities cause entrainment and water ingestion to commence very close to the vent source 24 . Our experiments consequently overestimate water entrainment and underestimate natural fountain heights. Nevertheless, our reconstructed MERs for D sw ≪ 1 in Fig. 5 are consistent with expectations on the basis of recent hydrovolcanic calculations by ref. 24, which conservatively limit water depths to be less than about 200 m for the inferred MER.
The particle volume fraction present in the jet mixture at the source is where V p is the volume of particles in the mixture and V tot is the total volume of the mixture. The particle volume fraction is used to calculate the bulk density of the mixture ρ 0 = (1 − ϕ 0 )ρ f + ϕ 0 ρ s where ρ f is the interstitial fluid density at the source and ρ s is the particle density. In addition to particle volume fraction, particle size (inertia) and, to a lesser extent, density determine the particle-fluid coupling regime that affects the gathering and sedimentation of particles from explosive eruption columns, ash clouds and related PDCs 19,22,[52][53][54] . Metrics for the particle-fluid coupling regime are the Stokes and sedimentation numbers: Here, τ p is particle response time to fluid accelerations, τ f is the characteristic fluid flow timescale and τ s is the particle settling time. Where St ≈ 1 and Σ ≈ 1, particle inertial and buoyancy effects modify the entrainment and sedimentation properties of eruption columns. CCF eruptions occur in collapsing regimes where the excitation of SWs at the top of the momentum-driven fountain region of the eruption column is governed by well-established mechanics of turbulent fountains 19,29 . The rise height of a multiphase fountain can, consequently, be predicted on dimensional grounds with the spatially averaged source momentum and buoyancy fluxes: where M 0 = πr 2 0 u 2 0 , B = πr 2 0 u 0 g ′ 0 and C is determined with the source parameters 19,[55][56][57] . Oscillations of the fountain about an average H ftn occur at a frequency: where C ftn ≈ 0. 5 (refs. 19,29). Both H ftn and f ftn are defined at the source with subscript 0 and at the air-water interface with subscript i. Assuming particles are not lost from the fountain mixture during rise through the water column 47 , the momentum and buoyancy fluxes of the fountain mixture entering the air layer can be estimated with measurements of the fountain radius and rise speed at the air-water interface.
The characteristic fountain frequency f ftn is distinct from the buoyancy frequency N, which indicates the response of the stabilizing air or water-column density stratification to inertial effects of the fountain 19,29 . Where f ftn > N, the frequencies of sediment waves descending along fountain margins will be determined by both f ftn and N, whereas where f ftn < N, the frequencies of these processes will be determined solely by the stratification, N (ref. 19).
For deep-water fountains, depending on jet-source conditions (Fig. 5), the relatively low frequencies and high amplitudes of fountain-top height oscillations are governed by f ftn and the frequency at which the ambient density stratification oscillates in response to the perturbation of the fountain, where ρ a,0 is the ambient density at the source height and the density stratification dρ a (z)/dz is assumed to be linear (Fig. 3a). We define N wtr for water layers, N air for air layers and N for combined water and air layers.
For shallow-water fountains breaching the water surface, there are three additional modes that can contribute fountain-height oscillations: (1) f air ftn defined for the subaerial part of the fountain; (2) the

Nature Geoscience
Article https://doi.org/10.1038/s41561-023-01160-z restoring stratification frequency N of the combined water and air layers, defined with the density difference between the water layer and air over the maximum fountain rise height through both layers; and (3) N air of the air layer defined using the density difference between the water layer and air over the subaerial fountain rise height (Fig. 3b). Low-frequency and high-amplitude fountain-height oscillations are governed by N wtr and N air whereas f air ftn governs relatively higher-frequency and lower-amplitude oscillations. Fountain-top oscillations are an expression of the build-up and collapse of periodic SWs around the upflowing fountain core (Fig. 3a,b). An additional mode is corner frequency related to inertial effects of the breaking waves and overturning motions characterizing stirring as well as entrainment and mixing across density interfaces as SWs descend: where r sw and u SW are the SW radius and speed as they descend next to the fountain. Usefully, we predict SWs to excite a distinctive and easily recognized fountain-top oscillation corresponding to this frequency (ref. 19). Fountain mixtures collapsing periodically from H ftn as SWs enter the water layer as momentum-driven jets with a characteristic jet entrance length scale 46 where Q sw and M sw are the volume and momentum fluxes of sediment waves calculated with their measured r sw and u SW as they enter the water layer. This jet entrance length scale will govern, in part, whether SWs are erosive or depositional on impact with the ground. The timing between SW impacts at the sea surface or seabed compared with the time for a SW to transform into a PDC and flow out of the SW impact zone will determine whether the spreading PDCs are fed intermittently or continuously. The ratio of an advective timescale for PDCs spreading out of the SW impact zone is τ PDC = r im /v PDC whereas the period between successive SWs is τ ftn = 1/f ftn . The ratio of these two timescales forms the 'collapse frequency' number In the special limit ℱ C ≫ 1, successive SWs coalesce such that periodic collapses typical of eruptions in collapsing regimes feed an approximately continuous mass flux to PDCs, resulting in massive or bedded deposits with negligible terracing. By contrast, as ℱ C → 0, periodic column collapses drive increasingly discrete PDC events. The MER of erupted material exiting the volcanic vent is a key source parameter used to compare eruption size and intensity and to initiate models of eruptions 25,36,58 , where A 0 = π[(r 0 ) 2 − (r 0 − Δr) 2 ] is the area of a ring vent geometry, with Δr the vent width, and is associated with caldera-forming eruptions 22 . Constraints on MER after an eruption are typically inferred from observations of eruption column height and estimates of total mass erupted divided by eruption duration 58,59 . Column height and MER estimated in this way, however, do not explicitly constrain unique values of -Ri 0 and ϕ 0 and, in turn, are not ideal for predicting whether eruption columns will collapse periodically or continuously. To link values of -Ri 0 and ϕ 0 associated with realistic natural MER values in Fig. 5, we use a mean caldera fissure vent radius r 0 = 4,000 m and fissure width Δr = 50 m, mean particle density ρ s = 1,500 kg m -3 and mean gas density ρ g = 0.2 kg m -3 and vary the source velocity over 100 ≤ u 0 ≤330 (refs. 21,35).

Data availability
Source parameter data for experiments conducted in this study are presented in Table 1

Code availability
Code used for spectral analysis of fountain-top height oscillations and deposit profiles are available upon email request: jgilchri@eoas.ubc.ca.