Experimental Multiblast Craters and Ejecta—Seismo‐Acoustics, Jet Characteristics, Craters, and Ejecta Deposits and Implications for Volcanic Explosions

Blasting experiments were performed that investigate multiple explosions that occur in quick succession in unconsolidated ground and their effects on host material and atmosphere. Such processes are known to occur during phreatomagmatic eruptions at various depths, lateral locations, and energies. The experiments follow a multi‐instrument approach in order to observe phenomena in the atmosphere and in the ground, and measure the respective energy partitioning. The experiments show significant coupling of atmospheric (acoustic)‐ and ground (seismic) signal over a large range of (scaled) distances (30–330 m, 1–10 m J−1/3). The distribution of ejected material strongly depends on the sequence of how the explosions occur. The overall crater sizes are in the expected range of a maximum size for many explosions and a minimum for one explosion at a given lateral location. As previous research showed before, peak atmospheric over‐pressure decays exponentially with scaled depth. An exponential decay rate of d¯0=6.47×10−4mJ−1/3 ${\bar{d}}_{0}=6.47\times 1{0}^{-4}\,\mathrm{m}\,{\mathrm{J}}^{-1/3}$ was measured. At a scaled explosion depth of 4 × 10−3 m J−1/3 ca. 1% of the blast energy is responsible for the formation of the atmospheric pressure pulse; at a more shallow scaled depth of 2.75 × 10−3 m J−1/3 this ratio lies at ca. 5.5%–7.5%. A first order consideration of seismic energy estimates the sum of radiated airborne and seismic energy to be up to 20% of blast energy. Finally, the transient cavity formation during a blast leads to an effectively reduced explosion depth that was determined. Depth reductions of up to 65% were measured.

An explosion-a sudden, rapid change of a material's volume that it imposes on its surroundings-forces that surrounding medium to rapidly compress such that the resulting pressure change does not propagate with the same speed as a smaller pressure change would. Small pressure changes can be described within the linear acoustic approximation, which assumes small pressure changes from a locally static pressure, and results in waves traveling with a (locally constant) speed of sound. Larger pressure changes cause adiabatic heating in air which locally increases the propagation speed and can lead to dramatic steepening of an initially smooth pressure wave into a discontinuity-a shock (Crighton & Scott, 1979;Garcés et al., 2013;Muhlestein et al., 2012). In an isentropic approximation (reversible process at constant entropy) a shock pulse has characteristic properties such as amplitude and duration that scale with the explosion's energy and the density of the medium in which the pulse travels (Kinney & Graham, 1985).
Scaling properties enable the establishment of phenomenological regimes that depend on scaled parameters, such as a scaled length. For example, for the depth d of a subsurface explosion, a scaled depth can be defined bȳ where E b is the blast's energy (Holsapple & Schmidt, 1980;Sonder et al., 2015). Using this method blasts of any energy may be categorized into deep, intermediate and shallow blasts. Deep blasts are contained in the ground and do not eject material (̄8 × 10 −3 m J −1∕3 ) . The host material's weight and strength are large enough to "contain" the blasts. Energy is dissipated by friction and anelastic alteration, or transported elastically as seismic waves. At intermediate scaled depths (̄≃ 4 × 10 −3 m J −1∕3 ) , material is excavated efficiently, which results in the largest craters. Shallow blasts (̄4 × 10 −3 m J −1∕3 ) create a smaller crater. They male smaller craters because more of the blast's energy couple with the atmosphere and fewer with the host, resulting in a large atmospheric pressure pulse. These regimes are backed up by extensive studies from military and mining research (Dillon, 1972;Ehrgott et al., 2011;Holsapple & Schmidt, 1980;Lee & Mazzola, 1989;Qiu et al., 2018), as well as research motivated by volcanology (Ambrosini et al., 2002;Goto et al., 2001;Ross et al., 2013;Sato & Taniguchi, 1997;Sonder et al., 2015;Valentine et al., 2012). For example, two blasts at same depth, which created craters with radii that differ by a factor 2, had energies that differed by a factor 2 3 = 8. Similar phenomenological regimes exist for a blast wave propagating in air. The distance from explosion source, r, may be scaled by blast energy and air density ρ̄= (2) The reference density ρ 0 is a value known from a case for which the scaled distance is known. Similar to ̄ , may be used to categorize an observation distance into far (̄6 × 10 −2 m J −1∕3 ) , in which the peak pressure drops with −1 , intermediate (̄≃ 6 × 10 −3 m J −1∕3 ) , or near ( 10 −3 m J −1∕3 ) (Kinney and Graham, 1985).
SONDER ET AL.

10.1029/2022JB023952
3 of 30 Less studied, from a volcanological perspective, are the effects of scaled depth on monitoring signals such as seismic, acoustic, and infrasound, particularly in cases involving multiple explosions occurring in rapid succession. Crater structures and ejecta products of such blasts are analyzed, and allow to connect their geometries and stratigraphy to energy, explosion locations and sequencing. These field findings also reveal the complexities of the natural processes, which limit the straight forward application of simple explosion models (Taddeucci et al., 2010). Some factors controlling the dynamic behavior and energy scaling have a common base with other applications of explosives in the fields of military or mining research (Ambrosini & Luccioni, 2006;Qiu et al., 2018). Such applications allow the scaling of lengths with a blast's energy, and use the depth below the surface to quantify its confinement. The scaling relationships were found experimentally, and while in detail the phenomena associated with a subsurface explosion depends on factors such as a host material's strength and density, rough phenomenological regimes can be identified that are primarily related to energy and depth combinations. Energy scaling was experimentally verified across length scales ranging from 10 −2 to 10 3 m, and energies from 10 3 to 10 15 J (Sato & Taniguchi, 1997;Strange et al., 1960;Vortman, 1968). Energies of most volcanic blasts fall into this range , and experimental research suggests that decompression driven blasting results in a similar scaled crater radius to depth relationship, even though changes in absolute values may have to be made, due to the different characteristics of the energy source (Liu et al., 2020). This motivates either direct applicability of the methods or a version adapted to volcanic activity.
Here we report results of experiments that focus on the effects of multiple explosions, closely spaced and timed, on ejecta, crater morphology, and geophysical signals. Such explosions show different behavior depending on the state of topography and host conditions at the time of explosion. Both conditions are varying rapidly, which causes ejecta jets to become asymmetric ( Figure 1 and video S1-S8 in Supporting Information S1), and can also be observed in the field on the volcanic scale (Voight, 1981). A volcanic explosive source was replaced by timeand energy-constrained chemical explosions. Previous experimental studies showed that this approach has important implications for field-scale analysis and interpretation (Bowman et al., 2014;Goto et al., 2001;Graettinger et al., 2014;Macorps et al., 2016;Sato & Taniguchi, 1997;Sonder et al., 2015;Valentine et al., 2014Valentine et al., , 2015. In these previous experiments explosive charges were detonated separately, and the effects of each single detonation on the surface morphology and ejected material were studied before detonating the next charge. Although that approach is relevant to many volcanic settings, observation shows that during explosive eruptions many explosions can occur closely spaced in time (Matoza et al., 2014;Park et al., 2021) or simultaneously, superposing their tephra jets, to create one single cumulative eruption column . Our study tests whether The example shows the jet of the third detonation in the "pad 1" configuration. Location markers for the third charge are bold. (See also videos S1 and S5 in Supporting Information S1).
SONDER ET AL.

Experiments
For each of the experiments reported here six charges were buried and detonated in test pads which were filled with unconsolidated granular material. The setup roughly follows previous studies on craters, each of which was created by more than one explosion ("multiblast craters") in which charges were detonated, and their blasts studied one at a time Sonder et al., 2015;Valentine et al., 2012). The use of unconsolidated granular material as host is motivated by the focus on the effects of multiple explosions. We note that any natural occurring rock material looses strength with scale (e.g., Hoek, 1999). Specifically in volcanic environments host material tends to be fractured. The presented experiments do not aim to determine the difference to blasting in solid rock material. The explosive material was Pentex™, which is a proprietary compound material with major components including trinitrotoluene (TNT) and pentaerythritol (PETN). It has a specific energy of 4.85 × 10 6 J kg −1 ; each charge had a mass of 90 g which corresponds to an energy of 4.37 × 10 5 J. The six charges were detonated in a timed sequence of 0.5 s between each detonation. Accuracy of detonation timing was better than 10 −3 s. This timing was selected to ensure that the ejecta jet of each blast interacted with that of the preceding blast. Two plan-view configurations were set up; one with three charge epicenters in a line ( Figure 2b); another with three epicenters corresponding to the apexes of a triangle (Figure 2c). Charges were arranged vertically on top of one another, at two depths, 30 and 60 cm. At the given blast energy 30 cm corresponds to a scaled explosion depth of 3.95 × 10 −3 m J −1/3 , a value very close to optimum excavation conditions. Horizontal spacing was chosen, such that the horizontal neighbor charge location would be within the footprint of a single blast at optimum depth, but close to its border. At pads 1 and 3 the upper charges were detonated in sequence, followed by the three lower charges. At pads 2 and 4, charges beneath each epicenter were detonated in a sequence of shallow-first and deeper-second ( Figure 3).
The blast sequences were monitored by high-speed and normal speed video cameras deployed on ground and on unoccupied areal vehicles (UAVs). A set of six cameras was arranged in a hemicycle, at a distance between 20 and 30 m to accurately capture directions of ejected materials. UAV-based video was recorded to determine lateral jet directions and material motion. High-speed cameras recorded at 300, 500 and 5,000 fps.
Seismo-acoustic records were made using a combination of seismometers, geophones, infrasound-microphones ("infrasound sensors") and higher frequency broadband microphones ("acoustic microphones"). The deployed seismometers and infrasound sensors fit into the SEED broadband category (band code "C", Ahern & Dost, 2012). Seismometers and infrasound sensors were recorded at 400 or 500 Hz. Deployed infrasound sensors had a flat frequency response between 3 × 10 −2 Hz and Nyquist frequency. Two types of the acoustic microphones were used, with linear (±2 dB) response from 3.15 Hz to 20 kHz and 4 Hz to 80 kHz (Table 1). Despite the familiar name "acoustic microphones" these sensors range far into the ultrasonic range. Recordings in this frequency range are rare for volcano seismo-acoustics, which often end around 5 kHz Nyquist . Some existing high-frequency field scale recordings are available from Goto et al. (2014) and Yokoo et al. (2002) who sampled at up to 50 kHz, and Ichihara et al. (2009), who report recordings at 200 kHz sampling rate.
From these sensors seismo-acoustic measurement stations were assembled for specific purposes. Station type A was dedicated to measure the radial decay of airborne-and ground-based blast signals. For each of the type A stations a 3-component seismometer, an infrasound microphone and two acoustic microphones were used. The seismometer was placed 1 m below the ground surface, and the infrasound sensor was just below the surface. The microphones were mounted 4 m above ground, pointing toward the blast source, and just above ground, pointing downwards. Seven type A stations were placed every 50 m in a radial line, starting at 30 m distance from the test pads center, so that the last station was at 380 m distance ( Figure 2). Station type B was dedicated to the depth dependency of blast signals. One station was assembled which consisted of three 3-component seismometers, placed 132, 75 and 18 cm below the surface, and one infrasound sensor, placed just below the surface. Station B had a distance of 30 m from the blast pads center (Figure 2). Station type C was dedicated to measure the angular dependency of the airborne signals. For each of them two acoustic microphones were placed 2.44 and 1.22 m above ground. Type C stations were placed in a 30 m radius semi-circle around the center of the blast SONDER ET AL.
10.1029/2022JB023952 5 of 30 pads. Angles range from 0° to 180° and were arranged so that the 90° station was also the start of the type A radial line ( Figure 2). Seismo-acoustic setup also included a line of 23 geophones to record ground speeds at 12-100 m distance along the type A radial line.
Ejected material was collected in two box arrays, separated at an angle >45° to collect material from 2.5 to 13.5 m from the charge assembly's center. The sample arrays were re-positioned for each experiment, so that they were always centered around an explosion site. One array was typically at an angle φ = 90°. The other array had different orientations for each pad, because other equipment and arrangements restricted the available space.
After the charges had detonated and ejecta jets had dissipated, photographs of the produced compound craters were taken for photogrammetry (structure-from-motion) analysis. Photographs were taken using (a) the same UAVs that also recorded blast videos, and (b) using a standard SLR camera, operated by a ground-based person. A subset of the photographs was the base for digital elevation models (DEMs) that were created using the commercial photogrammetry software Metashape™, generally following previous experiments . The resulting DEMs have a spatial resolution between 1 and 1.5 cm for pads 1-3, and 2.5 cm for pad 4. All crater profiles-and sizes presented below are based on these elevation models.

Seismo-Acoustics
An explosion is a sudden volume increase of a substance ("the explosive"), often caused by a highly exothermal chemical reaction. The expansion rate is larger than the speed of sound of the surrounding material which, in a general case, may be solid, liquid or a gas. The process compresses its environment faster than the latter can transport (at the speed of sound), which causes an oversteepening of the related pressure pulse and eventually produces a pressure discontinuity ("jump"). Close to the source the pressure jumps from ambient (atmospheric) value to a maximum and then relaxes back before sinking below ambient pressure and again relaxing back. At larger distances the propagation speed approaches the speed of sound and the pressure discontinuity relaxes to a steep, but finite slope.
Using features of wave-forms recorded by microphones and/or seismic sensors it is possible to estimate the blast's energy, provided that scaling laws assumed in such models are valid. The scaled peak pressure and scaled impulse of a blast in air depends on the scaled distance where the pressure is measured (Kinney & Graham, 1985). This relationship can be used to determine the scaled distance of each microphone record, and with that the energy of each blast wave can be estimated. This resource will be used as a reference model, and referred to as KG85 data (or -model). For these blasts in air, the main fundamental three quantities to be scaled are distance, time, and pressure. As in the case for underground blasts distances can be scaled with blast energy E b . Additionally, the relatively high atmospheric homogeneity allow further specification of the atmospheric density, which is often written in terms of transmission factors for scaled distance and time. Scaled distance, time, and pressure are given by Figure 2. Multi-sensor stations were placed in a radial line every 50 m starting at 30 m distance from the test pads. Each station included compact broadband seismic and infrasonic sensors as well as broadband ("acoustic") microphones. Acoustic microphones were placed in a 30 m radius semicircle around the center of the test pads. Another set of microphones was placed in a radial line from the test pads ranging from 30 to 80 m. 12 geophones were placed every 2.5 m starting at 12 m from the pads center, and 11 more along the same direction every 5 m following that. The last geophone had a distance of 99.5 m from the pads center. Six identical cameras recorded the experiments also in an arc of about 30 m distance. Other cameras recorded from a 50 m location. In pads 1 and 3 the upper charges (buried at 30 cm depth) were the first three to be fired, before the lower level (buried at 60 cm depth) was fired in the same lateral sequence as the upper ones. (c, d) In pads 2 and 4 charge pairs located at same horizontal location were fired consecutively (upper level before lower level). The main observation direction is shown by the large, gray arrow. The origin of the crater-specific calculations is set to pre-blast surface level for the vertical coordinate (z), and the geometric center of the charges (centroid) for the lateral coordinates (x, y indicated as dashed lines). Markers at the surface level show the "epicenters" above the charges. A circle shows epicenters above the first fired charge, a square shows an epicenter above the last fired charge, and a cross shows the epicenter above the in-between fired charges. SONDER ET AL. 10.1029/2022JB023952 where p a is the atmospheric pressure and the transmission factors f d , f t for distance and time, respectively, take the density into account in which the blast pulse propagates. They are given by Here, T is temperature and index 0 refers to values of a known blast case. The model only applies to explosive shocks in air. Energy E in Equation 3 refers to the blast's energy, E b , and assumes that the blast wave did not pass any material boundary. In case of subsurface blasts this value has to be reduced, since part of E b is dissipated and transported in the ground. If enough energy is not contained in the ground this part, E a , creates a (smaller) shock pulse in the atmosphere. Measurement of E a is described in Section 3.3.1.
Another widely used quantity to measure a blast's intensity, damage potential and energy is its impulse per cross-sectional area (Bush et al., 1946;Guzas & Earls, 2010;Kinney & Graham, 1985;Schnurr et al., 2020), which can be obtained as the time integral of the initial positive pressure peak of a microphone pressure curve as Here t s is the start time (time of arrival of the pulse at the sensor's location) and t 1 is the time of first zero crossing of the pressure curve ( Figure 6). This time interval always contains the peak pressure. The corresponding scaled impulse is a compound of scaled pressure and time components The KG85 data provides values up to a scaled distance of 3.1 m J −1/3 (500 m kg −1/3 ). According to this data set the scaled pressure and scaled impulse decay with 1∕̄ at relatively large distances ( 10 −1 m J −1∕3 , 20 m kg −1/3 ). The explicit values for the decay are. Each seismometer in any of the stations had three components, North ("N" or "1"), East ("E" or "2") and vertical ("Z"). Components were aligned vertically (positive downward, "Z"), radially (positive pointing away from the direction of the blast source, "N" , "1" → "R") and to the transverse direction (perpendicular to radial, "E" or "2" → "T"). b Vertical distance relative to local ground surface: positive above, negative below. Direction in parentheses is the direction of the microphone maximum sensitivity. c Upper limit refers either to Nyquist frequency or to sensor limit, see text. The ±1 dB frequency range of the 40 AE, 40 AO is 5 Hz-10 kHz; frequency range of the 40 BE is 10 Hz-40 kHz.
As is common in the analysis of blast waves (Garces, 2018;Kinney & Graham, 1985), peak pressures were not directly read as the maximum of the measured pressure curve, but impulse I 1 was calculated and compared to a function representing a blast pulse shape. We used a modified Friedlander shape see for example, Marchetti et al. (2013). Here p is the scaled peak pressure, and α is a positive number that determines the strength of the exponential's decay. The value of p p that fits the measured I 1 best was used for the peak overpressure. The impulse reference data are somewhat unclear, because the reference interpolation function (Equation B2) deviates from the reference data points by 17%. Since both, interpolation function and data points are part of the same publication (Kinney & Graham, 1985), we decided to follow their measured data instead of their interpolation parameter. Therefore one parameter of the interpolation function was adjusted, the constant ̄0 in Equation B2, which provides a better fit. Interpolation formulas, an explicit listing of all interpolation parameters and a graph that shows the difference between the two values of ̄0 are given in Appendix B, Table B1, and Figure B1, respectively.
Ford et al. (2014) determined distance-and depth-dependent energy partitioning of explosions above and below ground using a model for the airborne signal that, after some re-formulation (Appendix C), can be written as Here b 1 = 1.15 × 10 −7 s m J −2/3 and ̄3 = 1.2 × 10 −3 m J −1∕3 . Evaluated at ̄ = 0 this model expects a ca. 7% smaller scaled impulse (factor 2 −1/10 , ≃ 0.93) at a given distance compared to a free air blast. A larger discrepancy exists with respect to the KG85 data: The two constants for the −1 dependency, a I,ref , and b 1 differ by a factor 0.51. Determination of E a from the scaled radial to scaled impulse relationship, as done in Section 3.3.1 (Equation 15) using b 1 instead of a I,ref yields a factor ( ,ref∕ 1) 3∕2 ≃ 0.37 reduced values for E a . The data set presented here does not contain a zero depth or free air blast, and therefore cannot decide for one of the models. Energy values listed in Table 3 used the KG85 constant, and should be adjusted if used in connection with Equation 9.
Equation 9 and microphone records of previous blast sessions, carried out in very similar host materials and with similar explosives, show that scaled impulse decays rapidly with scaled depth (Appendix A). A somewhat more accurate match with experimental data is obtained for the peak pressure dependency on depth. Therefore the following is formulated using a peak pressure dependency. At depths ̄5 × 10 −3 m J −1∕3 peak pressure can be approximated by a product of an exponential which contains the depth part and an amplitude containing the radial dependency: Here the scaled depth related constant ̄0 = 5.4 × 10 −4 m J −1∕3 . This approximation is valid for scaled depths smaller than 5 × 10 −3 m J −1/3 ( Figure A1).

Energy Scaling and Crater Morphology
As described above, the scaled radius of a single sub-surface blast depends on the scaled depth of explosion. At several places, below analysis relies on a functional relationship of the two scaled values. Sonder et al. (2015) published an empirical function that uses understandable parameters for this dependency: The subscript 1 of this radius refers to the first blast in a sequence. The meaning of the four parameters is: ̄o pt is the "optimal" scaled depth of explosion; the scaled depth that crates the largest scaled crater of radius max .
9 of 30 0 is the scaled crater radius an explosion creates at 0 scaled depth. And b is a decay constant that determines how strong radius decreases with scaled depth at depths larger than ̄o pt . The values used for these parameters reflect the combination of unconsolidated host material and explosive charges, which were the same for Sonder et al. (2015) and the here presented experiments. Values used for Equation 11 are ̄o pt = 3.85 × 10 −3 m J −1∕3 , max = 7.51 × 10 −3 m J −1∕3 , 0 = 2.38 × 10 −3 m J −1∕3 and b = 1.72 × 10 −5 m 2 J −2/3 .
There is an upper limit for the size of a crater structure (volume or radius) which is created by a sequence of a (potentially) unlimited number of explosions at the same lateral location. Despite the (potentially) unlimited cumulative energy spent by the explosions, any single explosion can only break up finite amounts of material and eject material to a finite distance, thereby creating a "multi-blast crater" of finite size. The upper limit of a crater's size, for example, given by its radius, r ∞ , is the balance between amount of ejected material by an explosion and the amount of material falling back into the crater directly and collapsing from the crater walls. Quantifying the behavior from the first explosion in undisturbed ground to many explosions, Sonder et al. (2015) showed that the change of scaled crater size decreases with increase in the number of explosions where the independent variable n is the number of blasts that occurred, and n 0 = 0.9 is a constant, which was determined from experiments. n can be expressed as the ratio of total (cumulative) blast energy of a sequence to single blast energy n = E b,tot /E b , which is useful in cases of changing scaled depths or E b , so that n may assume non-integer values. For a sequence that starts from zero crater radius, Equation 12 has the solution where n 0 = 0.9 is an experimentally determined constant. We note that, because of the exponential asymptotic and n 0 < 1, approaches ∞ already for n = 3 or 4.
In subsurface blasting the depth of explosion (depth of burial) is usually measured as the distance to the horizontal surface. For volcanic applications most explosions occur in-or under a crater or under existing, non-flat topography, which removes the reference level from which depth was determined. Experimental investigation showed that in a first approximation, the shortest distance to the surface can be used as reference level (e.g., Valentine et al., 2012). For the case that an explosion occurs under the center of a crater, the crater bottom is the appropriate first-order approximation to measure the explosion depth. Taking into account the crater shape Sonder et al. (2015) showed that the mass distribution above the charge imposes a "confining force," and an "effective depth" can be defined from a crater's topography that reduces the scatter in crater size measurements.
The study also showed that the effective explosion depth rarely deviates by more than 10%-20% from the closest point to the surface.

General Observations
For all pads, the initial blast transported the greatest mass of material. From the main observation direction this charge was located at the top-left end of the linear setups of pads 1 and 2, and at the top-rear corner of the triangular setups of pads 3 and 4 ( Figure 3). Size and speed of these initial blasts (jets) had a range of 15-25 m, and speeds of several 10 m s −1 away from the ground. The initial ground motion exceeded 100 m s −1 . These values are comparable to previously conducted experiments (e.g., Valentine et al., 2012). Ejecta jets of the quieter blasts showed similar thinning behavior as was observed in previous experiments for blasts under pre-existing crater-topography Ross et al., 2013). Some jets had a main direction that was not vertical, but had a certain direction toward the main (temporally changing) crater void showing similarities with previously conducted off-center blast experiments . For pads 2 and 4, for which the lower charges were fired only 0.5 s after the upper charge (at same lateral location), the perceived loudness (not measured amplitude) of these lower charges was significantly larger compared to the previous optimum depth blast. In contrast, for pads 1 and 3, for which lower charges were fired 1.5 s after the upper charge at same lateral location, the blast noise was significantly muffled (Table 2).
SONDER ET AL.

Jets, Craters and Ejecta
Unlike past experiments in which a crater was analyzed after each individual blast, the timing of these multiblast experiments only allows for inspection of the final crater and ejecta. This crater is the cumulative product of six blasts that migrate vertically and laterally through the host. The blast sequences in pads 1 and 2 created craters elongated along the axis of the charges. The final craters of the triangular blast sequences (pads 3 and 4) were more round, with some visibility of single-charge crater outlines in the triangle's corners ( Figure 4).
The deepest points of the pad 1 and pad 2 craters were located between the central and right charge positions in the x-direction, and in close proximity to the symmetry line along the charges in y-direction. The lower right charge was always the last to detonate. The crater profiles preserved a stepped floor centered over the final charge ( Figure 4). The ejecta showed a prominent ray (ridge of material) that extended from the final charge location out of the crater in the direction of elongation (φ = 180°). Parts of the ray could be traced more than 10 m from the crater. For pad 1, one of the ejecta sample arrays was in line with this ray (video S1 in Supporting Information S1); in this direction the ejected mass per area was a factor ≃10 higher compared to the material collected by the array perpendicular to the charge line ( Figure 5). Also, mass distribution is better described by an exponential distribution in the φ = 180°-direction compared to the 90°-direction which is better approximated by a power Note. The "left" and "right" labels refer to the jet directions as seen from the main observation location. Polar-and inclination angles are also illustrated in Figures 2b and 2c. a Delay of the lower charges, relative to the upper charge at same lateral location (cf. Figure 3). b Polar angle is counted counter clock wise, and 0° along the axis parallel to the charge lines of pads 1 and 2, pointing to the right as seen from main observation direction. law. Isolated pieces of shallow-sourced gravel from pads 1 and 2 were observed further from the charges; one of them over 30 m away from pad 2, in the φ = 180°-direction.
The asymmetry of ejecta distribution around the linear charge array is similar to what was observed in previous off-center multiblast configurations with temporally well separated charge detonations . However, in those experiments a steep ejecta ring was formed on the side of the crater opposite to the direction of jet inclination . This steep ejecta rim was not observed in the here presented, overlapping blast sequences.
The triangular blast sequences of pads 3 and 4 produced more equant crater shapes resembling blurred circles around the triangular blast centers (Figure 4). Compared to the linear setups the deepest points of the craters were located laterally closer to the centroid and had a larger distance to the last blast's center. The pad 3 crater had a low point between the first and second (lateral) blast locations. Pad 4 had the low point close to its centroid. Both of the craters had shallow slopes in the vicinity of the crater rim (the topographic high around a crater center). Slopes were steeper closer to the center of craters. Ejecta were concentrated in three main directions for pad 3, and two for pad 4. Compared to the linear charge setups, the observed ejecta concentrations of the triangular Note. Only signals from the radial microphone line were used. All E a values were derived from a fit to the impulse-distance relationship (Equations 6 and 8). Only I 1 -values that followed an r −1 -dependency were used for the fit (Figure 8). For the loud blasts of pads 2 and 4 (blasts 2, 4, 6) the r −1 dependency ended for r > 130 m which was the case for microphones at distances up to 130 m ( ≤ 1.71 m J −1∕3 , 276 m kg −1/3 ). a Number of microphones used to fit the radial dependency to the data. b Minimum and maximum distance of the microphones used to determine E a . sequences were less pronounced. The ejecta concentrations originate from one vertex of the charge configuration to bisect the opposite side of the triangle (video S3 in Supporting Information S1). The pad 3 sequence had ejecta concentrations correlating to all three lateral charge positions. In the pad 4 sequence ejecta rays only correlated to blasts 3, 4 (φ ≃ 150°) and 5, 6 (φ ≃ 30°), since the first two blasts occurred in an effectively radially symmetric setting (blast 1 under flat topography, blast 2 under an approximately radially symmetric transient cavity).

Radial Dependency of Airborne Blast Pulse
Our recorded pressure pulses show most of the characteristic features of a free air explosion. This indicates that enough energy was not contained in the ground. It is thus possible to estimate the un-contained energy, E a , which created a shock pulse in the atmosphere. Comparison to the known yield of the detonation charges, E b , can then give information of the effect of explosion depth.
The more contained blasts did not create large enough blast pulses to make a reasonable comparison with the KG85 reference data. However, all initial and the perceived louder blasts of pads 2 and 4 (blasts 2, 4, and 6) created wave forms that were consistent with blast pulses and could be compared. In those cases peak pressure data were in agreement with a 1/r dependency at distances of up to 100 m. The impulse data stay consistent up to about 130 m distance (Figures 8a and 8b). At larger distances the values deviate significantly from 1/r.
To compare the measured impulse values to the scaled reference, an r −1 dependency was fitted to the un-scaled values of a given blast pulse, and the fitting constant a I was used to determine the location in the scaled graph. This determines an energy, E a ("atmospheric energy"), that creates the pressure pulse:̄ Since both distance and impulse scale with 1∕3 a , the procedure "moves" values on either axis when changing energy ( Figure 9). The result are scaled distances at the end of the KG85 reference scale ( 0.6 m J −1∕3 , 100 m kg −1/3 ). From the scaled distance a real distance r corresponds to the energy a = ( ∕̄) 3 , which can be interpreted as the energy not contained in the ground, and is smaller compared to the blast energy E b . E a was found to be around 1.5% of E b for the initial blasts, and about 5%-7.5% of E b for the perceived loud blasts in pads 2 and 4 (Figures 8c and 8d, Table 3).

Blast Energy, Charge Depth and Explosion Sequence
The first charge of a blast sequence detonated under a prepared, flat surface in unaltered host material. The following charges detonated under changed topography and somewhat altered host material. The lateral spacing (0.6 m, 8 × 10 −3 m J −1/3 ) corresponds approximately to the maximum crater radius for the blast energy, and similarly, the vertical spacing (0.3 m, 4 × 10 −3 m J −1/3 ) had, approximately, the optimum depth. Previous experiments showed that for such scaled distances the blast's jet changes shape and, if the topography above the charge has an overall orientation, it will also change direction (Ross et al., 2013;Valentine et al., 2015). We follow previous work (e.g., Valentine et al., 2012) and define the explosion depth, d, as the distance from the crater bottom to the charge. With this approximation, that is, neglecting the crater shape but not its depth, it is possible to evaluate Figure 5. Ejecta overview graphs for the four blast sequences ("pads"). Left column: Map view of layed out ejecta collection box arrays. Each small square represents a 0.5 × 0.5 m collection box. At larger distances parallel boxes were used to account for reduced ejecta load. The 2 m × 2 m square centered at x, y = 0 are the test pads, red crosses indicate the charge locations. φ-labels show the polar angle of the collection array. Large gray arrows indicate the main observation direction. Right column: plots of ejected mass per area at distances r from the crater center. The trends generally follow a power law, except for the 180° case of pad 1 and the 50° of pad 4, for which an exponential trend is a better fit. In those cases larger parts of the material in the main jet direction ended up in the collection boxes. Blue points show that there the decay power of the triangular pads (3 and 4) is slightly less strong than for the linear pads (1 and 2). Equation 10 for peak pressures of blasts that were shot at same lateral location for the two different blast delays (0.5 and 1.5 s) that were realized. For the pad 1 and 3 experiments this applies to the following pairs of blasts: (1, 4), (2, 5), and (3, 6). For the pad 2 and 4 experiments the blast pairs with same lateral location are (1, 2), (3, 4) and (5, 6). Evaluating Equation 10 for two peak pressures at same scaled distance leaves only the scaled depth to change. For example, considering the ratio of peak pressures of pad 2's blasts 2 and 1 relates the scaled depth of blast 2 to the previous one bȳ2 This formula can be applied to any of the above listed blast couples with consistent results (Figure 10a), showing that the so-derived depths are reduced by a factor 1.5-3, compared to their initial charge location relative to the surface. Since E b was the same for all blasts, the lower charge at the moment of its detonation can be estimated to be at a depth r =̄r b 1∕3 below the crater bottom at that time. And because the location of the lower charge is known to be 0.6 m below the original surface, the crater bottom can be estimated at z bottom = −0.6 m + d r (Figure 10b). The two delay times show that 0.5 s after detonation the crater bottom is deeper than at 1.5 s. At 1.5 s the crater bottom is about the same location that would be expected from a blast of energy E b at optimum depth.

Seismic Signal
We present here an initial estimate of seismic energy involved in the explosion experiments. A deep analysis of the seismic records will be part of future studies. The energy radiated from a radially symmetric seismic source may be estimated from the measured square velocity of the ground (particle) motion u r (e.g., Boatwright, 1980;Johnson & Aster, 2005) Here A and S are coefficients for signal attenuation and site response, respectively. ρ g is the ground density and c g the propagation speed of the ground, both at the observation location. For this first broad look at seismic energy these parameters are assumed to be constant. In this assumed energy estimate only one component of ground motion, radial component u r is non-zero. Other seismic components are therefore ignored in the following. Then E s can be approximated as In this approximation the proportionality factor F depends on a combination of ground properties and attenuation characteristics, but not on E s .
The multi-blast setting adds the difficulty that seismic signals originating from different blasts overlap at larger distances (e.g., for r ≳ 80 m, Figure 7c). From such distances only the cumulative seismic energy of a blast set can be determined: At closer ranges the blasts can be identified clearly in the u 2 signal. There u 2 decays quickly before the next pulse arrives, and integration over a finite time interval is a valid approximation for each blast (Figure 11a): When compared to the airborne signals, the seismic records show an inverted trend: The "muffled" blasts 1, 3 and 5 of pads 2 and 4, which had a much lower airborne signal created a larger seismic signal, when compared to  Comparison of peak pressure p p and impulse I 1 with respect to their applicability to estimate an explosion energy, and their compatibility to the scaled air blast data by Kinney & Graham, 1985 (KG85). (a, b) Lines show the best fit r −1 -dependencies for smaller distances that follow this trend. The impulse data show a better agreement with the r −1 -trend. Energies E a estimated from peak pressures are about a factor 10 smaller compared to the impulse-based estimates. The p p -values start to deviate significantly from the r −1 -trend at distances r > 100 m. blasts 2, 4 and 6 ( Figure 11a). This behavior serves as motivation for a potential energy partitioning scheme. For a given pad configuration the assumption is made that seismic and acoustic energy of a blast add up to a constant value. b = a + s + rem In this picture a change in E a of δE, for example, by a change of blast depth, would result in a change of E s by −δE. The remaining energy E rem stays constant. E rem may be viewed as a term that gathers the energy of all other parts of processes involved, including dissipation such as heat, friction and brittle deformation. The assumption that it is a constant means that changing the airborne energy, for example, by changing the blast's containment (depth), only changes seismic energy. The two "generalized elastic" terms can exchange energy, while the remainder term does not.
Energy conservation applies to each blast and to the cumulative case, which allow determination of the two unknowns F and E rem . With ΔE = E b − E rem the per-blast case becomes and the cumulative case is where 〈E a 〉 = ∑ E a,i /N b . The difference between the two cases is that for Equation 23 r is treated as independent variable, while in Equation 22 E a is independent. The average value 〈E a 〉 is a constant.
The left hand side values of Equation 23 were fitted to an r −2 dependency. The result shows the expected behavior: Pads 2 and 4 with the large airborne signals have smaller seismic signals when compared to their respective geometric counterparts pads 1 and 3 (Figure 11b). The per-blast data for the right-hand side of Equation 22 show a different trend of the pad 3 data compared to the other pads (Figure 11c). For small E a they are larger than the other pads, and then fall off quicker with rising E a . Since for the other pads no unique trend could be determined, pad 3 was treated separately, from pads 1, 2 and 4. For both cases intercept and slope were determined. Together with the cumulative case fit, values for ΔE and F were calculated. For pads 1, 2 and 4, ΔE was about 17% of E b (≃73 kJ), for pad 3 this value is about 10% (≃45 kJ). Highest values of E s are a factor two larger than highest values of E a . Consequentially in cases of observed higher E a blasts, seismic and airborne energies were comparable (Figure 11d). To be complete, values for F are 3.5 × 10 9 Js m −4 for pads 1, 2, 4, and 1.6 × 10 9 Js m −4 for pad 3.

Discussion
The sizes of the presented craters are larger than one-blast craters, but smaller than they could become when blasting many times (i.e., n ≳ 4) at these lateral locations with the same energy. Overlapping footprints from laterally shifting, time-separated explosions create compound craters with a footprint area that can be calculated from overlapping circles centered around blast locations . For a given explosion depth a radius is related to explosion energy by the scaled radius, and therefore the footprint area is, too. Evaluating Equation 13 for n = 1 relates the scaled crater radius after one blast r 1 to the scaled upper limit for many blasts ∞ . Solving for ∞ and multiplying by 1∕3 b yields the upper limit for this blast energy Here 1,max = 7.5 × 10 −3 m J −1∕3 is the maximum scaled radius of one explosion, which occurs at the optimum scaled depth and was measured in previous experiments . For repeated blasting at reduced depths an equivalent crater size, r ∞,red can be calculated by replacing r 1,max in Equation 24 with r 1,red . The measured Figure 11. Estimate of seismic energy from squared particle velocity. (a) Pad 2 test squared pressure signal of the infrasound sensor and squared particle velocity at first radial station (30 m distance). The seismic signal shows clearly identifiable pulses that can be separated into six time intervals. As described earlier for pad 2 the airborne pressure pulses of blasts 1, 3 and 5 are much weaker as those of blasts 2, 4 and 6. In contrast peak values of 2 are higher for blasts 1, 3 and 5, and somewhat weaker for blasts 2, 4 and 6. The trend is not as strong for the seismic signal as it is for the airborne signal. (b) Radial dependency of squared particle velocity integral. Measured values and fitted r −2 curves of Equation 23 are shown. Pads 1 and 3, with sequential shot depth configuration, produced a higher squared particle velocity integral, compared to pads 2 and 4 (interchanging blast depth). To a lesser degree, the triangular pads 3, and 4 had larger values when compared to the same shot depth configuration of the linear geometrical setups of pads 1, and 2. (c) Squared particle velocity integral dependency on E a . Despite some scatter, data from pads 1, 2 and 4 follow a common trend, while pad 3 data has a larger slope and offset. The black dashed line is a fit of Equation 22 to data of pads 1, 2, and 4 (correlation coefficient 0.67). The green dotted line to the pad 3 data (correlation coefficient 0.69). Cross markers show data form type A station closest to the blasts, circles data from the type B station. (d) Seismic energy plotted against acoustic energy for all pads. Black and green lines show the anticipated (linear) relationships using the derived values for F and ΔE. The second vertical axis shows E s relative to total blast energy E b . The elastic part is ca. 17% of E b for pads 1, 2, 4 and ca. 10% for pad 3. Gray dotted lines show the volcanic acoustic seismic ratio η = E a /E s (VASR, Johnson & Aster, 2005). The blasts had VASR values between 10 −2 and 1.
SONDER ET AL.
10.1029/2022JB023952 21 of 30 footprint radii fit into this picture: they range between 0.68 and 0.78 m, which is larger than the single explosion radius ( 1, max =̄1, max b 1∕3 = 0.57 m ) and smaller than the many-blasts limit. However, the crater sizes are not consistent with respect to the blasting sequence: in case of the linear setup, pad 1 (upper before lower charges) created a larger crater compared to pad 2 (interchanging charge depths), while in the case of the triangular setup pad 4 (interchanging depths) created the larger crater when compared to pad 3 (Table 4).
Equation 24 can also be used to estimate the final crater size of a hypothetical crater that would be the result of many blasts at reduced depth. It is then necessary to replace the maximum (scaled) crater radius with the reduced radius. The latter can be calculated from the scaled radius-scaled depth dependency for a single blast (̄ 1 ̄ r , Equation 11), using the scaled reduced depth value:̄∞ A footprint size estimated this way is larger than the measured two-blast crater, and ca. 7% smaller compared to the maximum possible crater (Table 4, Figure 12). Note. The reduced footprint is the maximum possible footprint when blasting at the reduced depth. The maximum footprint is the overall maximum that can be expected from the given blast energy Table 4 Measured-, Reduced-and Maximum-Expected Crater Sizes for the Tested Explosion Configurations 22 of 30 Goto et al. (2001) report that they conducted repeated blast experiments at same lateral location, and observed that the crater size only changed if the blast energy increased compared to the initial blast. This is an apparent difference to the results presented here and by for example, Sonder et al. (2015), Valentine et al. (2012), Ross et al. (2013), Graettinger et al. (2014) who report scaled sizes of multi-blast craters that are larger than that of a single blast. Goto et al. (2001) do not provide information on the explosion location, that is, whether the following charge was placed at same depth below the original surface, or at same depth under the crater bottom, or at some other depth. Therefore it may be permissible to speculate that the charges were placed at same absolute position (same distance to original surface) as the initial charge, which would correspond to an explosion depth reduced by the crater depth. Their reported data suggest that in this case the scaled depth would be lower than the optimal scaled depth that creates the largest crater, and therefore the repeated blasts could not increase the crater size.
Determination of the atmospheric energy E a from airborne impulse or peak pressure is possible for scaled distances up to about 5 m J −1/3 (800 m kg −1/3 ). At larger distances this type of analysis yields faulty values. A word of caution must be added, since the empirical models by Kinney and Graham (1985) and Ford et al. (2014) yield a factor 2 to 3 different energy estimates. A more in-depth analysis that focuses on the complete seismo-acoustic data set of the presented experiments may help here. For example, peak pressure of a weak shock (e.g., Muhlestein et al., 2012;Rogers, 1977;Young et al., 2015) decays with a power of radius slightly larger than 1. Such a dependency may be observed in the presented data ( Figure 8a). Other non-linear acoustic factors and near-field topography may also play a role (Maher et al., 2020). Nevertheless, both models evaluated here result in single digit values for the percentage of the energy ratio E a /E b . The relatively small amounts of explosives used have the advantage that the analysis does not have to deal with complications arising from drastically changing transmission factors (Equation 4), as in the case of large scale explosive events (e.g., Kim & Rodgers, 2016;Kim & Rodgers, 2017) or volcanic eruptions (Matoza et al., 2009).
For application to volcanic scale this validity range means that sensors would need to be deployed at a maximum range which depends on the explosion's energy and depth. Phreatomagmatic eruptions typically produce explosions with energies in a range of 10 9 J-10 13 J (240 kg TNT to 2.4 kt TNT, respectively) ). An explosion that occurs at a scaled depth of 4 × 10 −3 m J −1/3 , which creates the largest craters at a given energy and corresponds to non-scaled depths 2 and 86 m for above energy range, respectively. According to our results, at that depth about 1% of the energy would go into the airborne pulse, which corresponds to a maximum range of about 1-23 km, assuming 5 m J −1/3 as maximum applicable scaled range (Figure 8). For smaller events, sensors would have to be placed closer to a vent.
The changes in the apparent ("reduced") crater depth over time show that 0.5 s after detonation the crater is about a factor 1.5 deeper compared to 1.5 s after detonation. It is not clear whether this is the time of the transient cavity's maximum opening or not. The depth at 1.5 s is comparable to the depth of a single blast crater. For volcanic activity the timescale on which a crater forms is important. In this period part of the overlying mass confining magma in the ground is rapidly reduced, changing-or enabling non-steady state processes, such as magma-water mixing and phreatomagmatism (Büttner & Zimanowski, 1998;Lorenz, 1975) or decompression driven activity (Gonnermann & Manga, 2007). Assuming for a moment that crater formation duration scales with 1∕3 b , analogous to other blast related time and length (e.g., blast depth, crater radius), the presented results mean that for E b = 0.4365 MJ crater formation lasts on the order of 1 s, which corresponds to a scaled duration of 1.3 × 10 −2 s J −1/3 . An event creating a crater of about 15 m diameter would need 10 9 J  if created by a single blast, and would be formed in 1.3 × 10 −2 s J −1/3 × 10 3 J 1/3 = 13 s. A 25 m diameter crater would then need 44 s to form. It is, however, likely that other factors, such as variations in the host's material strength, complicate such a straight forward scaling approach. We note that previous work by Ohba et al. (2002) reported cube-root scaled "explosion cloud duration" dependency of scaled explosion depth of the same order. At optimum depth this scaled duration is 4 × 10 −3 sJ −1/3 , corresponding to ca. 0.3 s.
Despite such scaling difficulties the experiments show that hazards associated with explosions in the subsurface should be reevaluated. Multiple blasts may produce craters that allow deeper explosions to release energy at the surface, once they are no longer contained by excess load. Valentine et al. (2014) found that for well time-separated blasts, explosions at scaled depths equal or deeper than 8 × 10 −3 m J −1/3 do not erupt at the surface, but form subsidence pits, and have therefore a reduced hazard potential. The lower charges in this study were located 0.6 m below the original surface, which corresponds to a scaled depth of 8 × 10 −3 m J −1/3 , and would therefore stay contained in a single shot. The deeper lying charges of the pad 2 and 4 blast sequences not only erupted, but released significantly more energy into the atmosphere than the shallower charges, because explosion timing was such that the crater cavity was open. This effect corresponds to an effective reduction of scaled explosion depth by 65% 0.5 s after the upper charge and by 45% 1.5 s after the upper charge (Figure 10a). We note that such scaled depth reductions are equivalent to 22× and 5× increase in explosion energy (at same depth), respectively. The scenario of successively crater deepening, which is also of military interest (Antoun et al., 2003), cannot repeat indefinitely, since the following crater needs to move material from greater depth to the surface in a finite time window, which needs energy. More experiments are necessary to test where this limit lies, and what the exact crater formation duration is.
Analysis of the seismic signal reveals why the pad 3 crater is smaller compared to pad 4: Pad 3 had different attenuation-and coupling conditions leading to less energy available for seismic and acoustic pressure or momentum generation (ΔE), and more energy dissipated without momentum generation. The different coupling is the result of a variation in the pads host properties: On a subjective level, personnel preparing the pad for charge placement before blasting can confirm that pad 3 "felt" somewhat different compared to the others when punching holes for charge placement into the material. Such unintentional host variability highlights the sensitivity of the crater formation process to host properties (see also Macorps et al., 2016). The estimate of seismic energy and the energy partitioning analysis rely on good knowledge of E a . The assumptions made to estimate E s work well for large values of E a . At smaller E a (more contained blasts) scatter becomes larger, which suggests that the underlying assumption, that energy is partitioned only between seismic and airborne signal producing effects, does not apply there. The squared velocity and pressure signals of pads 1 and 3 emphasize this trend (Figures S9 in Supporting Information S1 and 11a). In a first order estimate the combined data from blasts in pads 1, 2 and 4 was between 10% and 20% of E b , and between 5% and 10% of E b for the blasts of pad 3. The experiments show how explosive energy is contained by friction, strength and inertia of the surrounding (overlying) material and how energy translates from driving ground-bound (seismic) to airborne processes, once the overarching containment parameter, scaled depth ̄ , is reduced.
Changes in host density and strength can make it harder to observe above phenomena on natural scale. Attenuation and site-specific response (Equation 17) depend on the host's strength and density, both of which also depend on the moisture content of the ground (e.g., Lamb et al., 1991;Rickman et al., 2011). Part of the energy which does not cause seismic or airborne waves is likely used to fragment the host rock (Dürig et al., 2021;Ouchterlony & Sanchidrián, 2019). Even though fragmentation is an important part of explosive volcanic processes (Zimanowski et al., 2003), the size of the presented experiments did not allow one to measure changes in collected representative samples of the ejected material. One-shot experiments in hard rock hosts show a strong decay of crater radius at scaled depths larger than "optimum" excavation depth (Rocco & Thomsen, 2010). If a blast sequence occurs in a hard rock environment the first blast likely fragments more material than following blasts (assuming similar energies), and consequently less energy would be available for other types of dissipation and elastic processes. Our experiments show that even in unconsolidated material large parts of a blast's energy is dissipated by fragmentation or friction.

Conclusions
Rapidly-timed subsurface blasts, occur in fields such as mining, geotechnical, military and medical applications (Arora et al., 2017;Mammadova et al., 2017;Qiu et al., 2018;Zhou et al., 2016). Our analysis of the ejecta, crater morphology, and seismo-acoustic signals should be applicable to those situations. We highlight volcanic eruptions, which commonly involve explosions in rapid succession (Dürig, Gudmundsson, Karmann, et al., 2015;Pistolesi et al., 2011). The results of this study provide insight on how to quantitatively interpret geophysical signals measured during such eruptions, as well as the resulting craters and deposits. The experiments show that energy is a robust parameter to relate the transient, dynamic phenomena, such as airborne and seismic pressure and stress waves and debris jets, with the long term products such as crater, subsurface deposits and ejecta. The observed changes of crater size due to the different charge sequence arrangements are on the same order as changes expected from different rock strength or water content (Nordyke, 1961;Rickman et al., 2011). Explosions from depths which, in well time separated settings, would stay contained in the ground, may breach the surface and erupt. They may breach if containment is temporarily reduced by a previous blast at smaller depth. The depth reductions observed in the presented experiments correspond to an energy increase up to a factor of 22. We emphasize that much of the presented physical signal analysis relies on (a) the high frequency records SONDER ET AL.

10.1029/2022JB023952
24 of 30 of airborne signal and (b) on the combination of relative near-field and far-field records. Deployment of such sensors hold promise for progress in seismo-acoustic volcano monitoring. Finally, the combined analysis of crater size and airborne signal showed the time dependency of the transient crater cavity. For the 0.44 MJ blast energy used here the transient crater is deeper 0.5 s after the blast than at 1.5 s. We think that this is important for the analysis of explosive volcanic blast dynamics.

Appendix A: Depth-and Distance Dependency of Peak Pressures From Previous Experiments
In previous blasting experiments Ross et al., 2013;Sonder et al., 2015;Valentine et al., 2015), a set of uncalibrated microphones was placed every 5 m starting at 5-30 m distance from the source. In all experiments the microphones were placed 10 cm above the ground facing toward the blast center. The blasts happened at various scaled depths with an emphasis roughly around optimum excavation conditions (̄≃ 4 × 10 −3 m J −1∕3 ) , but also deeper and some shallower blasts. Despite the uncalibrated pressure signal the raw signals were evaluated, since all sensors were of same model and therefore comparable. The result can be used to determine the relative depth dependency of impulse-and pressure signals, and compare them to other work (e.g., Ford et al., 2014). Signals were evaluated for peak pressure and impulse the same way as described for the here presented experiments in the main text.
Results show that the expected exponential depth dependency (Equation 9) underestimates both, pressure and impulse for deeper blasts ( Figure A1). Therefore a second term that only depends on scaled distance was added to the combined depth-and distance dependencies. At scaled depths smaller than 1.2̄opt ( ≃ 5 × 10 −3 m J −1∕3 ) the first term dominates, and the peak pressure show an exponential dependency ( Figure A1). At larger scaled depths peak pressures decay slower than this exponential predicts. More research is necessary, to clarify the slow decay. Bowman et al. (2014) suggest that ground motion dominates the airborne signal at larger depths. Best fitting values for the depth decay constant in the exponential is for the pressure case ̄0 = (5.4 ± 0.5) × 10 −4 m J −1∕3 , and for the impulse case ̄0 = (1.1 ± 0.3) × 10 −3 m J −1∕3 . ̄0 deviates by about 12% from the value found by Ford et al. (2014) responsible for depth decay (̄3 , Table C1).
We interpret this as good agreement for the range 0 ≤̄≤ 5 × 10 −3 m J −1∕3 . Figure B1. Effect of corrected value for ̄0 on the interpolation curve (Equation B2). For a reason not known to the authors the original value for ̄0 (orange curve) does not fit the KG85 data (black dots) well. We used a changed value, which better fits this data (blue curves, Table B1).

Appendix B: Interpolation Constants of KG85 Pressure and Impulse
The empirical equations for dependencies of blast overpressure, scaled impulse and scaled blast duration on scaled distance are as follows. Overpressure:̄=̄0 Values for the constants Z x,y are given in Table B1. For large distances, that is, the 1 in each of the factors in the above formulas becomes small when compared to the factor ∕ and can be neglected. Then and ̄ go with −1 :̄ 3.52 × 10 −7 s J −1/3 5.68 × 10 −5 s kg −1/3 Original value for ̄0 from Kinney and Graham (1985) is 6.61 × 10 −5 s kg −1/3 = 6.7 × 10 −2 bar ms kg −1/3 /1.01325 bar.  Appendix C: Impulse Depth-and Distance Dependency Ford et al. (2014) found the following model to fit scaled blast impulse, distance and depth: log 10̄ = 1 + log 10̄ + 3h − log 10 1 + 10 10 3h ∕10 (C1) Here h is the scaled height of burst, and energy was specified in kg TNT. Changing to scaled depth of explosion (̄= −h ) , this can be written as̄1 Constants b 1 and ̄3 are listed in Table C1. Ford et al. present the scaled impulse multiplied by ambient reference pressure, which is different from this study where scaled impulse is scaled overpressure integrated over energy-scaled time. We note that for ̄=h = 0 the depth dependent part reduces to 2 −0.1 ≃ 0.93, which is about 7% different from an exponential (e 0 = 1). For larger depths this difference is smaller, which justifies the use of an exponential depth part (Appendix A) without the reducing factor which is necessary above the surface:

Data Availability Statement
All measured data is hosted as datasets on Zenodo (Sonder et al., 2021). Data is also available separately as listed in the Supporting Information S1 document (data set S1).