PRESSURE WAVES AND TEPHRA DISPERSAL FROM VOLCANIC EXPLOSIONS : MODELS , OBSERVATIONS , AND INSTRUMENTATION by

Real-time study of erupting vents is important for both monitoring and scientific purposes; because direct in-situ study of erupting vents is impractical, our best tools for studying eruptions in real time involve monitoring eruptive products and waves that travel far from the volcano. The atmosphere is a particularly advantageous medium for studying propagation and transport of volcanic waves and products: acoustic waves pass through it with minimal scattering, particles follow predictable trajectories, and the atmospheric structure that affects both is well-monitored. Analyses of acoustic waves and tephra deposits can provide important information on eruptions including total explosive energy, volume, and fragmentation processes. Additionally, the hazards associated with these processes justify the need to understand and model them. Despite the apparent simplicity of volcanic-atmospheric phenomena, many open questions and difficulties remain. This dissertation aims to address some of the challenges and help develop a better understanding of volcanic-atmospheric phenomena. In this work, I discuss and demonstrate tools to improve our understanding of such phenomena. A general introduction to atmospheric physics and eruptive processes is provided in chapter 1. A particularly severe problem addressed by this dissertation is analysis of pressure waves from powerful volcanic explosions. Due to theoretical and numerical difficulties

associated with shock wave physics and the hazardous environment around exploding vents, existing theory, models, and observations are all insufficient to account for nonlinear shock wave propagation near the vent. This problem adds considerable uncertainty to potentially valuable acoustic inferences of eruptive activity. I address this problem in three ways. In chapter 2, I use numerical models of volcanic explosions to demonstrate a new framework for analyzing nonlinear pressure waves from powerful explosions, showing that tools developed for studying chemical and nuclear explosions can be adapted to study explosive volcanic eruptions. In chapter 3, I use existing acoustic theory and models to investigate an unusually powerful and wellinstrumented vulcanian eruption at Volcan Tungurahua (Ecuador), calculating the volume of erupted gas and tephra (∼ 0.5 km 3 ), classifying subsequent tremor into distinct mechanisms by its infrasound, and showing the relationship of volcanic lightning to vent activity. In chapter 4, I describe the development and use of a novel infrasound instrument (the Gem infrasound logger) intended to address limitations of existing instrumentation that particularly affect our ability to record shock waves. As the lowest-cost, lightest, and most flexible infrasound instrument currently available, the Gem is an ideal tool for recording shock waves in remote or hazardous settings where the risk of instrument loss must be tolerated and installation by drones with limited payload capacity may be necessary.
Finally, in chapter 5, I explore numerical modeling tephra transport from severe eruptions, focusing on two case studies. The first eruption, a 2015 lava fountain at Volcan Villarrica (Chile), produced a plume 6-8 km above the vent and deposited tephra in a narrow band extending tens of kilometers downwind. A custom Lagrangian model of tephra transport considering actual wind conditions at the time vii of the eruption shows good agreement with a map of the deposit obtained from field mapping and satellite imagery, including the finding of tephra grading perpendicular to the wind direction. In the second eruption, the 2013 vulcanian eruption at Tungurahua, I use the same Lagrangian model to calculate ballistic trajectories and times of impact with the ground, and show that coincident infrasound cannot be explained by other sources and probably originates in ballistic impacts. Infrasound due to ballistic impacts (which has previously not been documented) could be used to improve monitoring by enabling estimates of explosive properties to be made given ballistic properties, which could be detected and estimated rapidly after an eruption onset.  2.4 Peak overpressure ratio ( p−patm patm ) vs. scaled distance (z = rW −(1/3) ) for ideal shock waves produced by chemical and nuclear explosions of 1 kg T N T energy (from Kinney & Graham (1985)), figure 6-4. . . . . . 42 xvii 2.5 Arrival times of waves with a nonlinear propagation history relative to waves that have propagated only at the speed of sound (adapted from Kinney & Graham (1985)). Arrival time advance is calculated at long distance where the wave propagates at the speed of sound. Shock waves are only noticeably supersonic for a short period before decaying to approximately sonic speeds. . . . . . . . . . . . . . . . . . . . . . 43 2.6 Peak overpressure (P s = p−patm patm ) as a function of scaled distanceR = r(W/p atm ) −(1/3) and initial pressure for shock waves produced by compressed gas explosions (from Baker et al. (1978)). Note that the bar notation in the axis labels is different from the barred scales in sec- at Volcan Fuego (Guatemala), shown with synthetic traces found using numeric models of the volcano's topography and the inferred gas flux. Good agreement is seen between modeled and observed infrasound. Traces are normalized with respect to peak-to-peak amplitude of observed traces; vertical bars plotted with traces represent 2 Pa. Monitoring System (IMS) infrasound stations (Brown et al., 2014).
Gem power spectra as recorded (unscaled) and adjusted for decreased acoustic impedance at altitude (scaled) are presented (Bowman et al., 2017 in an eddy diffusivity of 750 m 2 /s). Point size and color again correspond to ballistic coefficient: small, blue points represent particles with small ballistic coefficients (generally finer tephra) and large, red points represent large ballistic coefficients (generally coarser tephra).
Due to diffusive effects of atmospheric turbulence, the deposit is noticeably wider than in the turbulence free model shown in figure 5.9.
However, similar effects of wind shear remain visible here: the axis of the deposit varies and ballistic coefficient grades from south to north. 162 5.11 Snapshots of tephra transport models at Villarrica. Left column assumes no diffusion (particles advected by steady winds only); right column assumes horizontal eddy diffusivity of 750 m 2 s −1 . Particles disperse more widely in the models including diffusion by turbulent eddies. However, tephra grading downwind to the southeast (due to large particles settling faster) and cross-wind to the northeast (due to wind shear) is evident in both models. The volcanoes of greatest interest to humans are found at the boundary between the solid earth and atmosphere, and many important eruptive processes are generated within the volcanic system and propagated by atmospheric phenomena. This dissertation explores these volcanic-atmospheric phenomena. The present chapter describes the eruptive and atmospheric processes that respectively create and propagate the volcanic phenomena studied in this dissertation-namely, pressure waves and tephra plumes. In chapter 2, I model shock waves radiated from volcano-like explosions and propose analytical methods for nonlinear infrasound. In chapter 3, I apply linear infrasound modeling and analysis to an unusually powerful and well-instrumented vulcanian eruption on 14 July 2013 at Volcan Tungurahua, Ecuador. In chapter 4, I discuss instrumentation challenges involved in recording pressure waves in the field and describe the development of a novel infrasound sensor-logger, which will facilitate recording pressure waves in the field in future studies. Chapter 5 describes the physics of tephra and ballistic dispersal in eruptions and the modeling of tephra fall in two very different eruptions (the 14 July 2013 eruption at Volcan Tungurahua and the 3 March 2015 eruption at Volcan Villarrica, Chile).

Background on explosion types
Explosive volcanism can occur in a variety of styles, all of which are driven by the behavior of pressurized, buoyant gas (the main species being H 2 O, CO 2 , and SO 2 ).
The processes that turn a magma's dissolved volatiles into accumulations of exsolved, pressurized gas that can explode at the magma surface depend on properties of the magma and volcanic system, which ultimately control the eruptive style.

Magmatic controls on explosivity
Volcanic gas begins as dissolved volatiles in magma; as magma approaches the surface, its pressure decreases, causing the solubility of volatiles in magma to decrease. Actual exsolution of volatiles depends on the availability of nucleation sites in addition to magma chemistry, pressure, and temperature (Proussevitch & Sahagian, 1998).
Once exsolved, gas behavior as bubbles depends on magma viscosity and bubble size: bubble coalescence and growth occurs more readily in less viscous magmas. Similar processes control bubble ascent: large bubbles in low-viscosity magma can rise buoyantly (with respect to the melt) with relatively little drag. Finally, the nature of the magma-air contact determines the vigor with which gas is released from the magma. A lava lake surface allows bubbles to burst immediately upon reaching the surface; conversely, a viscous or rigid seal prevents gas release until gas accumulates in sufficient volume and pressure to rupture it. Further, when the top surface of the magma is a downward-propagating fragmentation wave, relatively immobile gas bubbles burst continuously, leading to sustained gas and tephra emission (Alidibirov, 1994).
Coriolis parameter 2Ωsin(θ) (s −1 ) g acceleration from effective gravity (true gravity plus centrifugal force) (m · s −2 ) i the imaginary number mass flow rate of displaced atmosphere from a volumetric acoustic source (kg · s − Q any conserved quantity (any units) • degree of arc or latitude sr steradian An important characteristic of explosive eruptions is the fragmentation efficiency, the ability of the explosion to fragment magma into fine tephra. Like ordinary sediment, tephra is classified by size, the divisions being less than 2 mm (ash), 2-64 mm (lapilli), and greater than 64 mm (bombs if molten, blocks if already solidified) (Houghton et al., 1999).

Strombolian Eruptions
Strombolian eruptions occur in volcanic systems with low-viscosity magma that enables gas bubbles to grow, coalesce, and rise to the surface. Upon reaching the surface, the bubble pushes a membrane of lava upward, eventually bursting it and fragmenting it mainly into lapilli and bombs (Vergniolle & Mangan, 1999). Infrasound signals from strombolian eruptions are produced by bubble growth at the surface or gas release as the bubble bursts, and typically consist of a discrete pulse followed by a coda consisting of echoes (Witsil & Johnson, 2018) and crater resonance effects .

Lava Fountains
Lava fountain eruptions occur in conduit conditions similar to strombolian eruptions, but with a greater ratio of gas to melt. The resulting annular or dispersed flow transports clots of magma upward out of the conduit in a gas-dominated jet (Vergniolle & Mangan, 1999). Hawaiian activity can occur as "curtain of fire" eruptions along a fissure, or as a lava fountain erupting from a singe vent; the former style often transitions to the latter. Lava fountains have low fragmentation efficiency, and most tephra consists of bombs and lapilli that falls near the vent with ballistic trajectories, though convective plumes carrying finer material can also form above the lava fountains and carry material far downwind (Wolff & Sumner, 1999). Gas and tephra ejection is continuous in lava fountain eruptions, and so is the resulting infrasound .

Vulcanian Eruptions
When the top of a magma conduit is sealed by a dome or plug of viscous degassed magma, gas can accumulate until the seal ruptures catastrophically in a vulcanian eruption. Vulcanian eruptions produce compositionally and texturally diverse tephra, including blocks broken from the plug, material eroded from the conduit wall, and juvenile ash (Narvaez, 2014). Infrasound signals from vulcanian eruptions begin with a discrete pulse that can be quite powerful (e.g., Iguchi et al., 2008); the pulse can be followed by potentially long and complex continuous signals from gas and tephra emission. A gradual rise in pressure lasting around 1 s sometimes immediately precedes the main discrete pulse; this is interpreted to originate in pre-rupture plug deformation (Yokoo et al., 2009;Yokoo & Iguchi, 2010).

Plinian and Subplinian Eruptions
Plinian eruptions are characterized by a somewhat steady emission of gas and tephra from a downward-propagating fragmentation wave in an ascending magma column.
The fragmentation wave is a surface along which bubbles have grown enough to physically separate pieces of liquid magma, transitioning from a flow of melt containing bubbles to a flow of gas containing pieces of melt and tephra (Alidibirov, 1994). This transition removes the weight of the tephra from the magma column, causing the pressure to decrease, enabling further bubble growth and magma ascent.
The resulting sustained eruption can produce powerful and high columns that create extensive ash deposits; particularly powerful columns can penetrate the tropopause and disperse ash and aerosols into the stratosphere, where it can remain for months or years (Cioni et al., 1999). Infrasound from plinian eruptions is continuous and high-amplitude, corresponding to the continuous emission of gas and tephra (Matoza et al., 2011).

Eruptions Discussed in Subsequent Chapters
This dissertation focuses on two major volcanic eruptions. On 14 July 2013, Volcan Tungurahua, Ecuador experienced a powerful vulcanian eruption that produced pyroclastic density currents and was audible at least 180 km away (Instituto Geofísico-Escuela Politécnica Nacional (Ecuador), 2013). I analyze infrasound from this eruption in chapter 3 and model ballistic fall from the eruption in chapter 5. The second eruption discussed in the dissertation was a hawaiian eruption on 3 March 2015 at Volcan Villarrica, Chile, which produced a lava fountain reaching 1.5 km above the vent and a tephra plume that reached 6-8 km above the vent and drifted 400 km downwind (Sennert, 2015). I model tephra dispersal from this eruption in chapter 5. Additionally, I briefly mention a weak vulcanian eruption at Volcan Fuego (Guatemala) in chapter 4; infrasound recorded from this eruption is used to estimate eruptive gas flux and demonstrates the capabilities of a new infrasound logger.

Fluid dynamics of atmospheric air
Modeling volcanic-atmospheric processes requires an understanding of atmospheric gas dynamics. Classical mechanics and thermodynamics of air are governed by a system of partial differential equations enforcing conservation of mass, momentum, and energy (the Navier-Stokes system). It is common in fluid dynamics (and necessary, given the intractability of the full Navier-Stokes system) to neglect and approximate terms in the equations depending on the scale of the phenomenon being modeled. This approach enables us to model a variety of atmospheric processes with approximations that are highly accurate and easy to use. I use such approximations (described further in this section) to model pressure waves (chapter 2) and tephra (chapter 5).

Conservation Equations
The conservation of mass, momentum, and energy are perhaps the most fundamental tenets of classical mechanics (which includes geophysical phenomena like pressure waves and atmospheric transport). Local conservation laws are expressed in integral form as a balance between the rate of accumulation within a volume and the net flow out of the volume across its bounding surface: where Q is the conserved quantity per unit volume, φ is the flow vector of the conserved quantity per unit area, V is the volume being tracked, S is the bounding surface, andn is the outward-facing unit vector normal to the surface. When Q and φ are differentiable with respect to position and time, the divergence theorem can be used to put this law in the clearer differential form: The assumption of differentiability is not strictly satisfied in the presence of discontinuous shocks. However, for reasons of clarity I use this formulation for theoretical discussion, knowing that the more opaque integral form is required in the finitevolume numerical implementation (LeVeque, 2002). I use an Eulerian formulation of conservation equations (Chung, 2002) here and in chapter 2.
Importantly, conservation equations contain two unknowns: Q and φ. To close the system, a relationship must be established between conserved quantities and their fluxes. Mass flux can occur only by advection (φ = Qv), whereas momentum and energy flux can occur by advection, forces, or other processes. In chapter 2, the relevant fluxes are advection by particle motion (φ = Qv) and forces and work involving the pressure gradient ∇p. In chapter 5, the relevant terms are advection by steady winds that vary only with elevation (v = w(z)) and advection by turbulent eddies, which can be represented as diffusion (φ = −K∇Q), where K is a constant diffusivity tensor.

The Navier-Stokes System of Equations
The Navier-Stokes system, a set of three partial differential equations conserving mass, momentum, and energy, is the most general system of governing equations for fluids at macroscopic, non-relativistic scales. This system describes fluid phenomena at a broad range of scales, ranging from infinitesimal sound waves to planetary-scale circulation, and even fluid flow within stars and nebulae (Vallis, 2006).
The terms in the mass equation represent (in order) the change in density ρ balanced against the advection of mass (density times velocity vv).
The momentum equation includes (in order) the change in momentum (ρv) balanced against advection of momentum, the pressure gradient force, viscous forces τ , and body forces b like gravity, electromagnetic forces, centrifugal force, and the Coriolis force. ∂ρv ∂t Note that the momentum advection term ρv ⊗ v denotes ρ times the outer product (tensor product) of v with itself. In Cartesian coordinates, the tensor is expressed as Volumetric energy density is defined as E = E k + E i , the sum of kinetic energy density E k = 1 2 ρ|v| 2 and internal energy density E i = ρC V T (the product of density, specific heat at constant volume C V , and temperature). The energy equation includes the change in energy, advection of energy, work done against a pressure gradient, work done against viscous stress, work done against body forces, and conductive heat flow.
The Navier-Stokes system includes three equations but contains many more free variables. Therefore, additional equations are needed to close the system, including an equation of state relating pressure, temperature, density, and internal energy, and constitutive equations to calculate heat flow and viscous forces.

Body Forces in Atmospheric Fluid Dynamics
The main body forces relevant to Earth's atmosphere are gravity, centrifugal force, and the Coriolis force (Vallis, 2006). Although centrifugal force and the Coriolis force are "fictitious" forces, they appear in conservation equations constructed in rotating systems and are indistinguishable from "true" forces in that context.
The "true" gravitational force is the simplest of the body forces, equal to ρg grav , where g grav is a vector pointing toward the center of the earth with magnitude given by Newton's law of gravitation.
The centrifugal force and Coriolis force both depend on latitude θ and the rotation vector of the earth Ω (defined as having a direction pointing from the center of the earth to the North Pole, and a magnitude equal to the angular rotation rate of the earth |Ω| = Ω = 2π T (T being one sidereal day). Centrifugal force acts in the outward direction perpendicular to earth's rotation axis: F ce = ρΩ × (Ω × r e ) (where r e is the position vector relative to the center of the Earth). When working in a locally Cartesian coordinate system (x being East, y North, and z upward), the centrifugal acceleration can be added to gravity to make an "effective gravity" vector: g = g grav +Ω×(Ω×r e ). It is worth noting that even though the magnitude of the true gravitational force is much greater than the magnitude of the centrifugal force, g will generally not point exactly toward the center of the earth.
This does not cause problems on the spatial scales described in this dissertation: the vertical directionk is simply defined to be parallel to the effective gravity vector g so that the horizontal components of gravity are always zero.
The Coriolis force represents inertia perpendicular to the earth's rotation axis of bodies moving with respect to the earth and is defined as F co = −2Ω × vρ. This definition can be simplified by neglecting the vertical component of the Coriolis force, which is typically very small compared to effective gravity. In Cartesian coordinates, the horizontal Coriolis force is defined as F co = ρfk×v, where the Coriolis parameter f is defined as 2Ωsin(θ). The horizontal components of this vector can be written explicitly as (ρf v y , −ρf v x ). At large spatial scales in which the Coriolis force and the pressure gradient are the dominant horizontal terms in the momentum equation, circulation around low-pressure zones is anti-clockwise in the Northern hemisphere and clockwise in the Southern hemisphere.
Combining these terms, the body force vector may be defined as the sum of the effective gravity and Coriolis force b = ρ(g + fk × v). These terms are both relevant at long time scales only: they are negligible in typical volcano infrasound, but are essential in calculating mean winds that transport tephra (discussed in chapter 5) (Vallis, 2006).

Euler Equations
Although the Navier-Stokes system applies to fluids in general, this dissertation focuses on a single fluid: atmospheric air. Air has extremely small values of viscosity and conductivity; consequently, the governing equations can be simplified by dropping terms related to heat conduction and viscous stress (which do not affect waveforms at macroscopic scales, including when shocks are present). This leads to the Euler equations, a simplified set of conservation equations that assume inviscid, adiabatic flow.
The mass conservation equation is the same as in the Navier-Stokes system: change in density depends on mass advection.
The momentum equation is the same as in the Navier-Stokes system except that the viscosity term is neglected: In the energy equation, terms related to viscosity and heat conduction are neglected: To close the system of equations, an equation of state defining pressure in terms of ρ, v, and E is needed. Finding an equation of state can be facilitated by assuming the gas to be ideal. In an ideal gas, the share of volume occupied by molecules is negligible and molecules interact by elastic collisions only. The ideal gas assumption is most valid in temperatures that are not very cold and pressures that are not very high, and is a good approximation to atmospheric air at any realistic pressure/temperature condition at the surface of the earth and for the compressed gas explosions modeled in chapter 2. In an ideal gas, the pressure-temperature-density relation is p = ρRT (where the specific gas constant R is defined as R = C p − C V = (γ − 1)C V ). The ideal gas law allows E i to be conveniently rewritten in terms of pressure instead of temperature: E i = ρC V T = p γ−1 . (The specific heat ratio γ = C p /C V is an important descriptor of a gas's thermodynamics, and depends on gas chemistry: γ ≈ 1.667 in monatomic gases like He, γ ≈ 1.4 in diatomic gases the main constituents of air N 2 and O 2 , and γ ≈ 1.3 in many triatomic gases like the main volcanic volatiles H 2 O, CO 2 , and SO 2 ). The definition of energy suffices as a general equation of state (and is used in the shock wave modeling in chapter 2): In the specific case of isentropic disturbances (physical changes that do not change entropy), the simpler isentropic equation of state can be used: Eq. 1.11 is often referred to as the adiabatic equation of state, and its application to sound speed calculation results in what is often called the adiabatic sound speed. The term adiabatic is perhaps problematic, as it refers merely to processes that do not involve heat flow. However, eq. 1.11 requires the stronger assumption of a reversible adiabatic process that does not change entropy and is therefore called isentropic.
For example, rapid changes like shock waves increase entropy even in perfectly non-conductive gas, a situation that is adiabatic but not isentropic. Therefore, eq. 1.11 is valid over continuous parts of a pressure wave, but does not apply to discontinuous shocks themselves (Lighthill, 1978;Hugoniot, 1889).

Application of the Euler equations to pressure waves
The Euler equations without body forces can be justifiably applied to modeling pressure waves at macroscopic scales. Air has low viscosity and heat conductivity (Lighthill, 1978), meaning that viscous forces and heat conduction do not matter in ordinary atmospheric conditions at macroscopic scales. In addition to viscosity and heat flow, body forces can also be neglected over time and space scales relevant to pressure wave propagation, including when shocks are present.
In the immediate vicinity of a shock wave, viscosity and heat conduction do result in the transition between shocked and ambient states being finite but extremely thin (on the order of 10 −7 m) instead of the true discontinuity that would occur at the limit of zero viscosity and heat conductivity (Lighthill, 1978). The difference between a shock being infinitesimally thin instead of extremely thin does not affect shock propagation at macroscopic scales (LeVeque, 2002); further, because computational expense limits spatial grid spacing to values much greater than the shock thickness, the shock itself cannot be resolved in models. In particular, the precise values of viscosity, heat conduction, and shock thickness do not affect energy dissipation from the shock wave. Consequently, the Euler equations are used widely in theoretical work (Rankine, 1870;Hugoniot, 1889;Sachs, 1944;Taylor, 1950a) and numerical modeling of shocks in air (Bethe et al., 1944;Brode, 1959;Crepeau et al., 2001). Therefore, the numerical modeling in chapter 2 uses the Euler equations without body forces: (1.14)

Scale Invariance of the Euler Equations
Nondimensionalization is a powerful tool for understanding physical systems: dividing reference values from variables in governing equations can clarify how the system's behavior varies at different scales. The Euler equations can be nondimensionalized by substituting ρ =ρρ * , v =vv * , p =pp * , and E =ĒE * , and the derivatives ∇ =x −1 ∇ * and ∂ ∂t =t −1 ∂ ∂t * . In each substitution, starred quantities represent nondimensional variables and barred quantities represent characteristic scales of that variable; e.g.,ρ could be selected as the ambient atmospheric density, andx could be a length scale like the radius of a feature in the initial condition.
Nondimensionalization of the Euler equations results in the system ρ t which can be re-arranged as Barred quantities can be chosen arbitrarily, but careful selection of barred quantities can make each coefficient consisting of barred quantities reduce to one, simplifying the system. It is clear thatv (1.23) The formpρ is reminiscent of the sound speed c = γp ρ , suggesting thatv = c = p ρ ,p =Ē = γp 0 , andρ = ρ 0 are convenient reference values. The reference length x and reference timet are unconstrained except thatx =vt; any convenient length or time scale may be selected as long as their ratio is the velocity scale.
Having selected reference values, the nondimensional system is now scale-free: Consequently, solutions of the Euler equations at one set of scales can be adapted to different scales as long as the relationships among the reference values follow the constraints in eqs. 1.22-1.23. For example, the pressure waves found by modeling an explosion with initial radius r can be adapted to an explosion of initial radius 5r just by multiplying every possible length scale (e.g., sensor distance and wavelength) and time scale (e.g., wave arrival time and period) by 5 (applied in chapter 2). This is a technique is simple but powerful, considering that calculating solutions to the Euler equations is very computationally expensive.

Acoustic Waves
The Euler equations can be simplified immensely by assuming that pressure disturbances are much smaller than ambient pressure ((p 1 − p 0 ) p 0 ). Infrasound recordings are usually less than 1% of ambient pressure (Johnson & Ripepe, 2011), so this assumption is often justified at sites where sensors are placed but is generally invalid very near explosive sources (discussed in chapter 2. The simplification begins by taking the time derivative of the mass equation and the divergence of the momentum equation (1.15), and re-arranging: (1.28) Substitution leads to (1.29) Up to this point, the equation is still exactly consistent with eq. (1.15) in that no approximations or extra assumptions have been made. Now, the infinitesimal nature of the waves enables two approximations to simplify the equation.
First, without loss of generality, a reference frame velocity can be selected by Galilean invariance so that fluid velocity consists of infinitesimal perturbations from zero. Then, the final term in the equation is must be quadratic in velocity. Reasoning that products of two or more infinitesimal quantities are negligible, the final term can be dropped: This leaves a single equation of two variables, so another equation is needed to relate pressure and density and close the system. Because infinitesimal disturbances conserve entropy, the isentropic equation of state p/p 0 = (ρ/ρ 0 ) γ can be used. Differentiating the pressure-density relationship yields dρ = dp ρ 0 γp 0 (1.31) This is the linear acoustic wave equation, with sound speed c defined as c = γp 0 ρ 0 = γRT 0 (1.33) (Lighthill, 1978). The formulation in terms of temperature is found by applying the ideal gas law p = ρRT , where R is the specific gas constant (≈ 287 J kgK for air) and T is absolute temperature in Kelvins. Plane-wave solutions of the acoustic wave equation in a homogeneous space are of the form where A is the wave's amplitude, ω is angular frequency, and s = c −1î is the slowness vector (the vector oriented in the direction of propagationî with magnitude c −1 . Spherical waves take the form p(r, t) = A r e iω(t±r/c) . (1.35)

Geostrophic Balance
A different simplification of the Euler equations must be derived to explain large-scale atmospheric winds with periods longer than about an hour (for example, the ambient winds used to calculate tephra trajectories in chapter 5). At these large, gradually varying scales, the dominant horizontal forces are the Coriolis force ρfk × v and the pressure gradient, and the dominant vertical forces are gravity and the vertical pressure gradient (a situation referred to as hydrostatic balance: ∂p ∂z = gρ). These assumptions are justified at these scales except in the atmospheric boundary layer where interactions with earth's surface become relevant. Starting with equation 1.7, the Coriolis force is included as a body force, and terms including derivatives of non-pressure variables are dropped due to their gradual rate of change, resulting in ∇p − b = 0 (Vallis, 2006). From this balance, the east and north components of the wind vector can be written explicitly as The balance of the Coriolis force against the horizontal pressure gradient is referred to as geostrophic balance, and the resulting wind is called the geostrophic wind. Note that in geostrophic balance, the wind direction is perpendicular to the pressure gradient and therefore flows along isobars. Winds averaged over time scales greater than a few minutes (below which turbulence dominates; see figure 1.1) tend to be approximately geostrophic. With velocity defined explicitly using the momentum equation, the mass and energy equations can now be simplified. Geostrophic winds are much slower than the speed of sound (a regime in which air is relatively incompressible: ∇ · v = 0) and are perpendicular to the pressure gradient (so that v · ∇p = 0). Consequently, the mass and energy equations simplify to which are simply transport equations that describe the advection of a passive tracer.
So, in geostrophic flow, density and energy fluctuations (typically due to temperature or compositional variability) are simply advected with the wind and do not affect it (Vallis, 2006).

Earth's Atmospheric Structure and Processes
Earth's atmosphere is divided into five layers: the troposphere, stratosphere, mesosphere, thermosphere, and exosphere. Unlike the solid earth, the atmosphere lacks sharp material contrasts. Therefore, atmospheric layer boundaries (the tropopause, stratopause, mesopause, and thermopause) are defined by temperature maxima and minima, and therefore vary in elevation depending on season and latitude. Atmospheric flows originating on the surface of the earth usually remain within the troposphere due to the temperature minimum at the tropopause; only the most powerful plumes from volcanic or nuclear explosions reach the stratosphere as a result. However, acoustic waves from powerful explosions can propagate throughout the entire atmosphere and refract back to the surface at long distances from explosive sources .
Atmospheric properties including pressure, temperature, density, and wind velocity are central in modeling tephra dispersal (chapter 5) and pressure wave propagation (chapters 2, 3). The variability of these properties depends strongly on the time scale (figure 1.1). In particular, very little variation occurs in the band between a few hours and a couple tens of minutes, referred to as the "spectral gap" (Shuttleworth, 2012).
Therefore, it is convenient and common to represent atmospheric properties as the sum of long-period ("mean" or "advective") and short-period ("turbulent" or "eddy") components, where the latter is considered unresolvable and handled statistically.
Through most of the atmosphere, the mean wind flow is geostrophic (horizontal and controlled by the balance of the horizontal pressure gradient and the Coriolis force; see section 1.2.6). Accordingly, these winds flow mainly horizontally and parallel to isobars. The gradual variability of these winds (and of other atmospheric

Turbulence
Atmospheric turbulence mostly arises from air flow over an uneven topography and is destroyed by forming successively smaller eddies that are eventually dissipated by viscosity upon reaching very small spatial scales (Shuttleworth, 2012). Turbulent motions appear random and can be described as having a characteristic time scale τ ≈ 1 min ( fig. 1.1), a variable velocity variance σ 2 v , and zero mean velocity. Random motion of this sort causes properties advected with the flow to spread out, away from regions of high concentration toward regions of low concentration, similar to molecular diffusion. In fact, the diffusion equation describes this process well: where K is a tensor representing diffusivity by turbulence (instead of by true molecular diffusion) and C is the conserved atmospheric quantity being advected. When K is isotropic and constant over space, the diffusion equation simplifies to ∂Q ∂t = K∇ 2 Q. When modeling tracers affected by turbulence in a Lagrangian framework, the equation of motion consists of a decaying inertial term and random motions: where N (0, 1) represents a normally distributed random number with mean zero and variance one.

Inversion of Discrete Volcano Infrasound
Infrasound from volcanic explosions provides scientists with important information on explosion intensity, which is most easily estimated as the volume of air displaced by the eruption (the subject of chapter 3). Early work (Johnson et al., 2003) on this topic used a simple acoustic model including a linear monopole source and spherical spreading: q being the mass flow rate of displaced atmosphere, r the distance between vent and receiver, c the atmospheric speed of sound, p the infrasound excess pressure, and Ω atm the solid angle subtended by the atmosphere around the vent (2π for a half-space, and between 2π and 4π for a conical volcano). Importantly, this model neglects effects like topographic scattering, refraction by atmospheric winds and sound speed heterogeneity, and source anisotropy.
Development of a fast, user-friendly finite-difference acoustic model (Kim & Lees, 2014) permitted effects of topographic scattering and atmospheric refraction to be predicted easily. This advance enabled subsequent work to deconvolve these effects from recordings and improve accuracy of source volume estimates Fee et al., 2017;Kim et al., 2015). Recent work has even constrained multipole source components by recording infrasound on tethered balloons anchored over volcanic vents (Jolly et al., 2017).
Propagation effects discussed so far are purely linear-that is, they can be expressed exactly as a function to be convolved with a source function to obtain the recorded infrasound: G(t) being the Green's function representing propagation effects. Obtaining the source function is then simply a matter of deconvolving the Green's function from the recorded infrasound. When infrasound is recorded simultaneously at multiple sensors in different locations, a linear system can be set up and inverted to find the ordinary least-squares solution: where H is the generalized inverse of the matrix of Green's functions Kim et al., 2015).

Analysis of Continuous Volcano Infrasound
In addition to explosions, volcano infrasound can arise from continuous vent processes such as resonance of open-vent systems and jetting. Resonance occurs when a vent and conduit radiate particular frequencies in response to a source; it can be excited by both continuous (e.g., roiling of a lava lake surface) and discrete (e.g., strombolian explosions) disturbances. Vent resonance has been used to estimate vent geometry (Richardson et al., 2014;Goto & Johnson, 2011), explain tremor (Fee et al., 2010b) strombolian explosion coda (Witsil & Johnson, 2018), and infer lava lake stage .
Volcanic jets are rapid, sustained, turbulent emissions of gas and tephra from a vent. Although turbulent acoustic sources are difficult to model from first princi-ples, the problem of turbulent jet acoustics has received considerable empirical study due to the importance of noise reduction from jet engines. Volcanic jets produce infrasound whose spectrum can be related to properties of the jet like its diameter and speed . Jet noise spectra are often complex but can be fit approximately by models corresponding to large-scale turbulence or fine-scale turbulence. Additionally, supersonic gas jetting creates a type of tremor called "crackle", which consists of asymmetric pulses consisting of strong compressions followed by gradual rarefactions; crackle has been observed at Nabro volcano  and possibly Tungurahua , discussed in chapter 3).

Slowness vector and source location estimation using infrasound arrays
An additional application of volcano infrasound is to locate the absolute position or the direction of acoustic sources, particular when it is necessary to distinguish between multiple volcanic vents (Ichihara & Matsumoto, 2017) or track moving sources like pyroclastic flows (Ripepe et al., 2010), rockfalls (Johnson & Ronan, 2015), or lahars .
The ability of a set of infrasound sensors to locate sources depends on its geometry.
Two or three-dimensional arrangements of microphones installed much closer to each other than to possible sources are called arrays and can be used to identify the slowness vector of incident sound waves. In a two-dimensional array, only the two slowness vector components in the plane of the array can be calculated directly; the third component must be found by assuming a sound speed c and using the constraint Powerful volcanic eruptions produce shock waves when high-pressure volcanic gas is released into ambient air; the shock waves propagate nonlinearly as they exit the conduit and diffract out of the volcanic crater, and eventually decay into linear infrasound waves. The physical processes that occur in volcanic shock waves are poorly understood and much more difficult to model than ordinary infrasound. This problem limits the field of volcanic infrasound because in order to infer source processes from infrasound recordings, all propagation phenomena must be accounted for-including shock wave processes.
This chapter addresses the problem of shock wave propagation from volcanic explosions by combining analytical methods used to understand chemical and nuclear explosions with numerical modeling. Simple but powerful analyses have been developed to infer sources of such explosions; the most powerful tools include scaling analyses by which properties of recorded pressure signals can be related to source energy. However, such explosions occur at much higher pressure (and therefore, a much higher initial volumetric energy density) than volcanic explosions, and lack the severe topography around which volcanic shocks propagate; therefore, the methods developed for chemical and nuclear explosions cannot be applied directly to volcanic explosions.
To assess the utility of similar scaling analyses in the context of volcanoes, I ran numerical models simulating volcano-like explosions at many reasonable combinations of pressure and temperature, both with and without topography, and analyzed results of the models in the context of scaling laws. Results show that scaling laws are valid for explosions with no topography and are sometimes valid when topography is present. A primary control on the validity of scaling laws with topography is wave diffraction around the crater rim: short waves diffract less readily around the crater and are therefore recorded with spuriously weak amplitude along the ground surface.
Although it would be much more difficult than recording along a volcano's slope, recording infrasound either above or within the crater could avoid this problem.
Numerical models allow a wide range of initial conditions to be tested, but are time-consuming to set up and computationally expensive; scaling laws vastly decrease the number of models that need to be run. This combination of several pre-generated explosion models with scaling laws allows explosive energy to be estimated from an infrasound waveform. Although the numerical models used here are simplified compared to true eruptions (in particular, they do not account for the gas species and solid content of the erupted material), this provides a promising framework for infrasound analysis of discrete eruptions.

Introduction
As discussed in chapter 1, infrasound recordings can provide important constraints on source physics, including the important parameter of explosive volume. However, any conserved quantity (any units)  Assuming wave propagation to be linear is justified when pressure amplitudes are much lower than atmospheric pressure, which is normally valid at the relatively distant sites (hundreds of m to tens of km) where infrasound is typically recorded.
However, a pressure wave that is linear at the site where it is recorded may have a partly nonlinear history. Explosive source pressures can be as high as several tens of atmospheres (Clarke et al., 2015), meaning that pressure waves near cannot be assumed to propagate linearly near the vent in explosive eruptions. The potential existence of nonlinear propagation effects complicates infrasound analysis: they are difficult to model, cannot be implemented as a convolution, and cannot be removed by solving a linear system. Therefore, they are typically neglected, even in significant vulcanian eruptions.
Past research has provided essential theory for understanding the formation and structure of shock waves when propagation is nonlinear; however, we still lack simple analytical solutions analogous to equation 1.42. Further work has helped mature our understanding of shock waves in the special case of synthetic explosions when simplifying assumptions can be made (Kinney & Graham, 1985). However, a proper understanding of shock waves in the context of volcanic explosions still eludes us, despite recent attempts to study them with numerical models (Morrissey & Chouet, 1997), field data , and laboratory experiments (Medici et al., 2014).

Nonlinear Effects on Waves
Nonlinear effects on waveforms generally fall into three categories. First, areas along the wave profile with high pressure normally have a higher sound speed and particle velocity than other parts of the wave, meaning that they tend to propagate faster and overrun the rest of the wave. Consequently, the front of the wave tends to steepen into an extremely narrow shock, a process referred to as the wave oversteepening, breaking, or shocking-up. At high pressures, this process can occur even when the initial condition has gradual spatial variations only. Finally, the narrow, high-gradient shock front causes irreversible compression, meaning that the wave continuously dissipates energy by increasing the entropy of shocked gas.

Waveform steepening in nonlinear waves
Acoustic waves propagate at the isentropic sound speed plus the particle velocity.
When wave amplitude is not infinitesimal, these quantities vary along the waveform: compressions will have faster particle velocities and higher sound speed than rarefactions.
The method of characteristics (Lighthill, 1978) graphically illustrates the effects of variable propagation speed in a simplified one-dimensional scenario (although the outcome holds qualitatively in higher dimensions too). For particle velocity v and sound speed c = γp ρ (where p is pressure, ρ is density, and γ is the specific heat ratio introduced in section 1.2.3), the method of characteristics finds that certain quantities (referred to as Riemann invariants) remain constant from the initial condition along lines in x−t space (referred to as "characteristics") with slopes of v−c, v, and (plotted in fig. 2.1) v + c. Of those characteristics, v − c and v + c correspond to advected pressure waves (which can carry differences in pressure and particle velocity normal to the shock), while the characteristic with slope v corresponds to pure advection (which can only carry shock-parallel velocity and temperature-related density differences).
When v and c vary spatially, characteristics diverge in some regions and converge in others, causing the waveform to change in shape: compressions tend to develop sharp onsets and gradual drop-offs (fig 2.1).
However, the method of characteristics as so far implemented breaks down when characteristics intersect, and predicts unphysical multi-valued solutions (in the real world and in FV models, a location can have only one pressure). Therefore, a dissipative shock surface is needed to truncate characteristics to prevent them from crossing.
In order to truncate both the ambient and shocked-state characteristics, the shock speed must be c 0 < c s < (c 1 + v 1 ), where the unshocked air has sound speed c 0 and velocity v 0 = 0, the shocked air has sound speed c 1 and velocity v 1 , and the shock moves at c s . By invoking the Rankine-Hugoniot theory of jump conditions (discussed in section 2.2.1), one can identify a unique speed at which the shock must propagate (equation 2.7) as well as a means of dissipation-irreversible compression-particular to the shock (demonstrated in figure 2.3).

Rankine-Hugoniot Shock Jump Conditions
I have so far described the governing equations in terms of modeling wave propagation through the atmosphere. It is also sometimes necessary to explicitly describe the shock itself in terms of the shock's speed and fluid density, pressure, and velocity upstream and downstream of the shock. The problem of theoretically predicting changes across a shock in an inviscid, non-conducting gas was solved by Rankine   A: pressure waveform of passing shock. B: pressure-volume plot, color coded to match waveform segments in A. Black lines are isentropes. Initially, the air is at ambient pressure on an isentrope corresponding to low entropy (blue dot in B, first blue line in A). The arrival of the shock wave increases entropy as it irreversibly compresses the air (following the Rankine-Hugoniot relations), driving the state to high pressure on a higher-entropy isentrope (red). Subsequent changes are isentropic, and the subsequent rarefaction (green) and re-equilibration (blue) follow the new adiabat. After the shock passage, the air is again at ambient pressure but with lower density, higher temperature, and higher entropy than before.
(1870) (deriving jump conditions that conserved mass, momentum, and energy but increased entropy) and Hugoniot (1889) (deriving the non-isentropic pressure-density relationship in instantaneous changes). These works established the key physical principles of shock waves: that mass, momentum, and energy must be conserved everywhere, and that the ordinary isentropic equation of state p/ρ γ = constant is true everywhere except at the shock (meaning that entropy increases as the shock passes). Because the shock wave problem lay at the intersection of two esoteric, lateto-mature fields (thermodynamics, and calculus of discontinuous functions), it was solved slowly and controversially, and remained widely misunderstood even decades after the contributions of Rankine and Hugoniot (Salas, 2007).
Calculation of conditions across a shock can be simplified by adopting a reference frame moving at the same speed as the shock. In this reference frame, ambientpressure gas (represented with the subscript 0 ) flows at supersonic speed into the shock. At the shock front, the pressure, density, and sound speed of the gas increase and the gas velocity decreases, making the velocity of the shocked gas (represented by the subscript 1 ) subsonic. Conservation of mass, momentum, and energy following the Euler equations requires the following (referred to as the Rankine-Hugoniot relations) to hold: where the mass density E = 1 2 ρv 2 + p γ−1 is the same as in eq. 1.9, and velocity v is assumed to be perpendicular to the shock. These equations can be applied directly to standing shocks, with supersonic gas leaving a nozzle being a common example.
However, when studying explosions, the ambient air is typically stationary and the shock is not, so switching to a reference frame with stationary ambient air is necessary.
Further, it is convenient to solve for the jump conditions as a function of pressure (an easily and commonly measured quantity). With these ends in mind, the following relationships can be derived from the stationary-shock equations (with the subscripts 0 and 1 still referring to ambient and shocked states): The speed of the shock front is Note that c 0 < c s < (c 1 + v 1 ) and that c s is not a sound speed, fluid velocity, or their sum.
Notably, the Rankine-Hugoniot relations show that the shock speed, fluid velocity, and pressure ratio across the shock can be arbitrarily high; however, density ratio across the shock has a strict upper bound: ρ 1 ρ 0 < γ+1 γ−1 . The decoupling of density from pressure makes the shock's compression irreversible so that entropy necessarily increases across a shock; consequently, the final temperature of ambient air after reequilibrating to ambient pressure is necessarily higher than the initial temperature before the shock's arrival.
In atmospheric air, we normally assume that normal dissipative processes (viscosity and heat conduction) are negligible at macroscopic scales (due to air's very low values of viscosity and thermal conductivity). However, these dissipative terms are proportional to gradients of velocity and temperature (eqs. 1.3-1.5), which become extreme in the nearly-discontinuous shock. Consequently, these otherwise negligible processes control the width of the shock (approximately 10 −7 m in air). However, macroscopic features and processes of the shock wave are insensitive to changes in viscosity or conductivity as long as their values remain very small (Lighthill, 1978): in fact, the Rankine-Hugoniot jump conditions (eq. 2.5-2.7) apply when viscosity and conductivity are arbitrarily close to zero. In practice, numerical methods used to model shock waves-for example, the Lax-Wendroff-Leveque method employed by Clawpack (LeVeque, 2002)-typically include small diffusive terms for the purpose of numerical stability (Chung, 2002). Consequently, viscosity and heat conduction do not need to be modeled explicitly in theoretical (Taylor, 1950a) or numerical models of shocks (Bethe et al., 1944) except when investigating the very small spatial scales at which those properties are relevant (< 1µm).

Shock wave theory for chemical and nuclear explosions
A great deal of past research has explored shock waves in the context of deliberate and accidental explosions due to their importance in industry and warfare (Kinney & Graham, 1985). The resulting body of shock wave theory has helped elucidate shock wave propagation, particularly when shock waves initiate in simple and compact sources.
The most important aspect of an explosion is the total released energy W (referred to as the "yield"); waves from explosive sources that are equal in yield but different in other ways (e.g., pressure, temperature, volume) will resemble each other at sufficiently long distances. This makes an observed or modeled blast wave for an explosion of a given energy broadly applicable to other explosions of that energy. W can be measured in J or (more commonly with chemical and nuclear explosions) with kg T N T (defined here as 1 kg T N T = 4.184 M J, though this definition varies among sources).

High explosives
High explosives are defined as chemicals that decompose in a chemical reaction that propagates faster than the material's sound speed by mechanical energy transfer, referred to as a detonation. By contrast, the chemical reaction in low explosives moves through the explosive material at a subsonic rate and transfers energy by thermal processes; this process is referred to as a deflagration (Kinney & Graham, 1985). A high explosive and low explosive of equal energy might produce blast waves of equal intensity at a great distance, but the high explosive will be more effective at breaking an adjacent object (a property termed brisance) because it releases energy instantaneously before explosive material can be dispersed by the explosion. Additionally, high explosives are easier to model because it is sufficient to assume that energy release occurs homogeneously within the original explosive volume, rather than having to model the chemical reaction and explosion kinematics simultaneously.
High explosives, exemplified by TNT, have been the subject of considerable numerical and empirical study due to their industrial and military importance and the relative simplicity of their explosions (Kinney & Graham, 1985). Among the results of this work is a reference blast wave profile showing the evolution of a blast wave emanating from a free-air explosion of 1 kg of TNT in typical sea-level atmospheric conditions ( fig. 2.4).
Following the Rankine-Hugoniot relations (eq. 2.5-2.7), shock overpressure corresponds to shock propagation speed: shocks with higher overpressures propagate faster. This means that a powerful shock wave could arrive at a receiver in noticeably less time than a sound wave from the same source. However, because shock waves decay so rapidly, most of that arrival advance occurs from high wave speeds very near the source, and explosions must be quite powerful for supersonic propagation to noticeably affect arrival times ( fig. 2.5). Like other time scales involved in the Euler equations, arrival times can be scaled following eq. 2.9.
Although the high explosive case is more compact than low explosives or compressed gas explosions, the mass and finite radius of the explosive are significant very near the source. This differs from nuclear explosions, which are approximately point-like in that the explosive mass and radius are insignificant considering the immense yield. However, nuclear explosions are themselves complicated by the fact that considerable energy goes into electromagnetic radiation rather than the blast wave, and that air's equation of state changes at high temperature due to ionization and dissociation of molecules (Kinney & Graham, 1985).

Compressed gas explosions
The term "compressed gas explosion" generally refers to sudden releases of gas stored at pressures well above one atmosphere but below initial pressures of high explosions.    Kinney & Graham (1985)). Arrival time advance is calculated at long distance where the wave propagates at the speed of sound. Shock waves are only noticeably supersonic for a short period before decaying to approximately sonic speeds.
relatively little study compared to high explosives (whose applications and hazards are more widespread). Compressed gas explosions typically have much lower initial pressures (and therefore energy density) than high explosives, meaning that they are less point-like and have lower energy density. Because of effects of finite explosion radius, blast waves from explosions of equal energy will differ from high explosives near the source but tend to converge at long distance (figure 2.6). Consequently, every initial gas pressure must be modeled separately, though scaling laws (below) still apply to explosions of equal pressure and different volumes (Baker et al., 1978;Esparza & Baker, 1977;Blanc et al., 2017).

Scaling laws for explosions
Simple scaling laws allow a blast wave profile for one explosion to be scaled up or down to explosions of similar type but different energy. Scaling laws can be used to rescale simulations of explosions to new initial conditions that differ only in scale.
This criterion requires that topography or atmospheric structure either be absent or that its geometry scale along with the explosion. As shown in section 1.2.4, the Euler equations are exactly scale-invariant as long as velocity scaling (typically sound speed) is equal to length scaling divided by time scaling. Therefore, shock waves produced by an explosion of a sphere of compressed gas should scale up exactly if the only parameter that changed was the radius of the gas. The additional scaling tools discussed below provide approximate scaling of shock wave length scales with the cube root of explosive energy (W = (p exp −p atm )V exp /(γ−1) for compressed gas explosions).
Because the cube root of the explosive gas volume V exp is proportional to the length scale, scaling by energy is equivalent to scaling by length (and is therefore exact) for explosions of equal pressure, but is only approximate when explosive pressure differs Figure 2.6: Peak overpressure (P s = p−patm patm ) as a function of scaled dis-tanceR = r(W/p atm ) −(1/3) and initial pressure for shock waves produced by compressed gas explosions (from Baker et al. (1978)). Note that the bar notation in the axis labels is different from the barred scales in section 1.2.4. Each profile corresponds to a different initial pressure; the profiles converge at long distances. among explosions.
The scaling parameter k, related to the cube root of explosive energy, is central to explosion scaling. Three schemes with different degrees of generality can be used to calculate k. The first assumes equal atmospheric pressure and temperature among explosions: k = W (1/3) , commonly expressed in units of (kg T N T ) 1/3 (Hopkinson, 1915;Cranz, 1926). The second allows variable atmospheric pressure but assumes constant temperature: k = (W/p atm ) (1/3) (Sachs, 1944), which can be expressed in units of length (e.g., m). The third, which can be expressed in units of m −(1/3) s −(2/3) allows variable atmospheric pressure and temperature: k = (W/ρ atm ) (1/3) , where density ρ atm can be obtained from pressure and temperature following the ideal gas law p = ρRT (Kinney & Graham, 1985). Any scheme is valid as long as it is used consistently and its assumptions are satisfied, though the second scheme may be preferred for dimensional reasons (it has simple units and can be used directly as a length scale) assuming that temperature does not vary much.
To demonstrate the application of scaling laws, consider two explosions (with subscripts A and B ) differing only by their scaling parameters k A and k B , and a "reference explosion" with k ref = 1. Pressure-related properties of explosions, like the peak overpressure of a shock, can be defined for all explosions with similar but scaled initial conditions as a function of scaled distance r/k only: Consequently, if the peak overpressure profile of explosion A p A (r) is known for all distances r, then the peak overpressure from explosion B can be predicted at any given distance r B : Shock wave parameters with dimensions of time, like arrival time and positive phase duration, are only slightly more complicated: for some time-dimension shock Then, shock wave scaling among explosions can be determined as Obtaining large numbers of blast wave profiles is difficult and inconvenient, whether by experiment or numerical model. However, with the scaling relationships, a result from a single explosion can be generalized with simple arithmetic to other explosions of widely varying energies, making the scaling relationships a powerful tool in blast wave study of chemical and nuclear explosions (Kinney & Graham, 1985); however, their utility has not yet been assessed in the case of volcanic explosions, which are complicated by topography and comparatively low, variable explosive pressure.

Strong shock approximation
"Strong shocks" are waves whose pressure amplitudes are much greater than ambient pressure (Taylor, 1950a), enabling additional simplifications of their behavior. Strong shock waves in a homogeneous atmosphere obey the following proportionalities: In these equations, r is the radial distance traveled by the wave, t is the time since the explosion, c s is the shock speed, p 1 and p 0 are the shocked and ambient pressure, v 1 is the shocked particle velocity, and ρ 1 and ρ 0 are the shocked and ambient gas density.
Strong shock theory assumes that the wave is far enough from the source that effects of finite explosion volume are small but near enough that shock pressure is much greater than ambient pressure. Specifically, the shock pressure must be high enough that the density equation in the Rankine-Hugoniot relations (eq. 2.5) can be approximated by its infinite shock strength limit ρ 1 ρ 0 = γ+1 γ−1 . These assumptions are satisfied over long distances for nuclear explosions; consequently, blast wave recordings from nuclear explosions agree well with strong shock theory (Taylor, 1950b).
Strong shock assumptions are satisfied over a much smaller range for high explosives (Taylor, 1950a) and not at all for typical compressed gas explosions (including volcanic eruptions). However, they provide an illuminating upper bound on the rate of decay of a shock's overpressure and speed.

Volcanic Shock Wave Studies
Most previous shock wave research has studied chemical, nuclear, and (to a much lesser degree) compressed gas explosions. This past research, discussed in the preceding sections, has developed a body of theory and observations that make shock wave study possible despite the intractable governing equations. Such explosions differ substantially from explosive volcanic eruptions, and assessing the applicability of such tools to volcanic settings is a focus of this chapter. However, past work has worked on volcanic shocks specifically with modeling, lab experiments, and field observations. Nonlinear propagation should occur whenever a wave's pressure amplitude is nonnegligible with respect to ambient pressure. Vulcanian eruptions exemplify this condition: source pressures in vulcanian eruptions can reach 10-15 MPa (though are usually less than 5 MPa) (Clarke et al., 2015), compared to ambient pressures of 0.1M P a or less (depending on vent elevation). Nonlinear propagation must occur in these eruptions, at least very near the vent. Despite the certainty of near-vent nonlinear propagation in powerful eruptions, the phenomenon's effects being negligible has often been assumed or stated (e.g., Buckingham & Garces, 1996;Kim et al., 2015;Fee et al., 2017).
A few studies have attempted to elucidate nonlinear volcanic infrasound. Morrissey & Chouet (1997) first modeled volcanic shocks with the objective of recovering burst pressures. Their sophisticated model varied the particle concentration and burst pressure in the exploding gas using a system of eight partial differential equations (the Navier-Stokes system with two phases). However, aspects of the model were unrealistic: it neglected topography, assumed the source volume (normally unknown) to be fixed, and treated the atmosphere as consisting of steam rather than air. Despite these limitations, their results show a clear dependence of shock strength and shock speed on particle concentration and initial gas pressure. A more recent nonlinear acoustics project has applied the Lattice Boltzmann method-a Lagrangian formulation of subsonic nonlinear gas dynamics-to volcano infrasound (Brogi et al., 2018).
This work modeled the wavefields that would result from a pulse of subsonic gas into the atmosphere. It found good agreement with linear acoustic theory when sources are relatively weak (consistent with the assumptions of linear theory) but that sources with higher Mach numbers of injected gas (M a > 0.3) could create anisotropic radiation patterns. Further, it found that asymmetric waveforms like the Friedlander waveform, which are sometimes interpreted as originating in nonlinear processes (e.g., Marchetti et al., 2013) Visible pressure waves ("flashing arcs") have been observed in many eruptions (Perret, 1912;Nairn, 1976;Genco et al., 2014). Visible waves are sometimes attributed to condensation of atmospheric water vapor due to compression in the shock ; however, water vapor actually condenses in powerful rarefaction waves-not compressions-which trail the shock and move no faster than the speed of sound (Yokoo & Ishihara, 2007). The compressive shock front (which may be supersonic) may manifest itself as a change in the air's index of refraction (Perret, 1912) or as vaporization of pre-existing clouds (Yokoo & Ishihara, 2007).
Observation of supersonic propagation of flashing arcs would be valuable because it would enable recording of near-vent volcanic shocks from a safe distance. Unfortu-nately, many studies of flashing arcs either fail to measure the wave speed, or find a wave speed indistinguishable from sound. As an exception, Ishihara (1985) did record videos of supersonic flashing arcs at Sakurajima volcano, Japan, but did not record infrasound for comparison. More recently, Marchetti et al. (2013) observed features analogous to flashing arcs in thermal videos of strombolian eruptions at Yasur volcano, Vanuatu. This study tracked a cold feature in their thermal images apparently moving at supersonic speeds, though with high uncertainty in the feature's speed.
The interpretation of the wave being supersonic is problematic because the estimated wave speeds in different eruptions did not vary with infrasound amplitudes recorded farther away (waves that were more powerful at a distance should have propagated faster near the vent). Further, the waves did not decelerate as they spread outward from the vent (shock speed depends on overpressure following eq. 2.7, so decaying shock waves necessarily decelerate). Additionally, when the 1/r spreading rule is used to estimate wave amplitudes 20-80 m from the vent (where the thermal feature was tracked) from the relatively low pressure amplitudes recorded on distant microphones, the estimated overpressure of the wave is insufficient for it to be noticeably supersonic.
Although the pressure waves in these relatively small explosions were not significantly supersonic, the ability to observe pressure waves near the vent with thermal video is promising for future studies of more powerful eruptions.
The problem of volcanic shock wave propagation has also been addressed with lab models. Controlled shock waves can be generated using devices called shock tubes that contain regions of compressed and ambient gas separated by a diaphragm.
Breaking the diaphragm causes three wave-like features to form: a compressive shock surface that propagates at supersonic speed into the ambient gas; a contact disconti-nuity surface that moves at the fluid velocity and separates regions of equal pressure and velocity but different density; and a rarefaction region that propagates at variable but subsonic speeds into the compressed gas. The scaling properties of the Euler equations make small shock tubes a viable means of conveniently studying volcanoscale shocks inside a lab. Medici et al. (2014) used shock tube experiments to simulate a sudden release of pressurized gas from a conduit into the free atmosphere. Médici & Waite (2016) expanded on that work by using shock tube experiments to explain the occurrence of multiple shocks in a single explosive event. Simple explosions were already known to commonly produce waveforms containing either one shock (e.g., the Friedlander waveform) or two shocks (e.g., the "N wave", which has a trailing shock formed by overexpansion, collapse, and rebound of exploding gases in addition to the primary shock); however, they generally do not form complex series of subsequent shocks or pulses. Médici & Waite (2016) found that supersonic jets trailing the main shock wave in a shock tube experiment could form subsequent pressure pulses, analogous to pressure pulses emanating from volcanic jets. Chojnicki et al. (2006) used shock tubes for a different purpose: understanding the effects of suspended ash on shock propagation. Volcanic eruptions often release two-phase flows including volcanic gas and ash; dynamics of ash particles differ strongly from gas and complicate the physics of flow and wave propagation. This paper showed the inadequacy of theoretical predictions of two-phase flow in experiments analogous to vulcanian eruptions, and proposed a new set of empirical laws for those conditions. Methods I use numerical modeling to explore shock wave behavior in a volcanic context.
Unlike linear acoustic modeling discussed in chapter 3, nonlinear pressure wave models must consider the full Euler equations instead of the much simpler acoustic wave equation, and are therefore more complex, computationally expensive, and prone to instability. Consequently, I use a set of numerical methods, implemented in the code Clawpack ("Conservation LAWs PACKage"), to model shock waves from explosions of compressed gas (LeVeque, 2002). The scenarios I simulate are analogous to discrete volcanic eruptions, but have been simplified to facilitate modeling. As volcanic shock modeling progresses, future work will use fewer simplifications and model volcanic eruptions more faithfully.

Finite Volume Models
Clawpack uses the Finite Volume (FV) method, a class of numerical schemes used to solve partial differential equations. The FV method is similar to the Finite Difference (FD) and Finite Element (FE) methods, but has certain properties that make it wellsuited for modeling hyperbolic systems of conservation laws like the Euler equations.
FV models implement systems of conservation laws, normally expressed in the differential form ∂Q ∂t where Q is volumetric density of a conserved quantity (such as mass, momentum, energy, or charge) and φ = f (Q) is the flux of that conserved quantity (having dimensions of the conserved quantity divided by area and time). However, as mentioned in chapter 1, the differential form is actually found from the infinitesimal volume limit of the conservation law's integral form where V is a control volume with bounding surface S. When spatial derivatives φ are defined everywhere, the divergence theorem can be used to to rewrite this integral as However, this assumption is not met when shocks are present, so it is necessary to use the integral form of eq. 2.17 in which the conserved quantity is tracked over a volume instead of at a point.
The integral formulation of conservation laws lends itself well to the FV numerical scheme, which tracks a large number of control volumes (termed grid cells) that completely fill the model region; the conserved quantity within each volume is calculated in each time step according to the fluxes, and fluxes are calculated on the surfaces between cells according to the conserved quantity within each cell. As a result, any flux across a boundary is subtracted from the upstream cell and added to the downstream cell, meaning that the total conserved quantity within the model region never changes-a desirable numerical property when modeling conservation laws.
(Appropriate exceptions to this rule occur when waves encounter nonreflecting model boundaries.) For simplicity, I briefly discuss the discrete formulation of the FV method in one dimension. Higher dimensions follow identical principles but have more cumbersome notation. Time and space must first be discretized: the current moment in time is labeled n and is separated from the subsequent time step by ∆t, and the cell of interest is labeled i and is of length ∆x. The change in notation from t to n∆t and x to i∆x, with i and n being integers, reflects the discretization that has occurred.
In one dimension, Q n i represents the average density of a conserved quantity per unit length within cell i at time n (note that superscripts represent the time index and not exponents). So, the total conserved quantity within a cell is ∆xQ n i and the change within the cell between times n and n + 1 is written as ∆x(Q n+1 i − Q n i ).
Because fluxes are calculated between cells and not within cells, the spatial indices are somewhat different: the flux from cell i − 1 to i at time n is written as F n i−1/2 . Perhaps inconsistently, the time index n (instead of n + 1/2) is retained even though the flux occurs between times n and n+1; this is because the flux is purely a function of Q n and is not at all influenced by the later state Q n+1 . The FV method can therefore be concisely written in one dimension as ∆x(Q n+1 This numerical method generalizes to multiple spatial dimensions and multiple conserved quantities; although the notation becomes increasingly complex, it consistently maintains the principle that the change in a conserved quantity within a cell equals the net influx of the quantity into the cell from other cells.

Riemann Solvers
Notably, I have not yet included any information from the Euler equations except that they are conservation equations, and so the fluxes F remain undefined. The next step is therefore to define what a given flux F n i−1/2 must be given Q n . Unfortunately, the answer given by a simple average of the flux function of the two cells-

Non-Hyperbolic Terms
So far, I have discussed conservation laws without considering source terms. A more general formulation of a conservation law takes the form ∂Q ∂t +∇·φ) = S, where S(x, t) is a source term. In FV modeling, the concept of a source term is often interpreted liberally to encompass any non-hyperbolic term, including terms representing geometric spreading, second-derivative terms like viscosity (e.g., Alghamdi et al., 2011), external forces, and actual injections of new mass, momentum, or energy into the model. In the models discussed in this chapter, the only source term used is that of geometric spreading: the models represent propagation into three dimensions, but exploit symmetries in the problem to simulate the propagation using only one computational dimension (when spherical symmetry applies, such as a free-air explosion) or two computational dimensions (when cylindrical symmetry applies, such as an idealized conical volcano).
Geometric source terms arise from the definition of divergence in spherical or cylindrical coordinates. In the case of spherical symmetry (where ∂ ∂ϕ = ∂ ∂θ = 0) being modeled in one computational direction (ther direction, so that φ = φ rr ), the divergence can be rewritten as ∇ · φ = 1 r 2 ∂ ∂r (r 2 φ r ) = ∂φr ∂r + 2 r φ r . In the case of cylindrical symmetry in cylindrical coordinates (where φ = φ rr +φ zẑ ), the divergence is rewritten as ∇ · φ = 1 r ∂ ∂r (rφ r ) + ∂φz ∂z = ∂φr ∂r + ∂φz ∂z + 1 r φ r . By contrast, the divergence in one dimension is merely ∇ · φ = ∂φr ∂r , and in two cartesian dimensions is ∇ · φ = ∂φr ∂r + ∂φz ∂z . Therefore, when a single computational dimension is used to represent either two or three dimensions, ∇ · φ will consist of the normal ∂ ∂r φ r term present in the one-dimensional divergence as well as an extra term 1 r φ r or 2 r φ r . This nonhyperbolic term representing geometric spreading can be considered a source term.

Accuracy of Numerical Methods
All discrete numerical methods are approximations to continuous methods, and their accuracy is often described in terms of the order of the approximation. For example, a simple linear approximation would consist of a constant plus a term of degree one; such an approximation would be referred to as first-order and its errors would be of second order. A second-order approximation would include a constant and terms of degree one and two, and its errors would be of third order. The order of an approximation determines the nature of the errors it results in. Odd-order approximations suffer from dissipative errors in which the signal is artificially smoothed, resulting in energy loss (e.g., fig. 2.7a-b). Even-order approximations suffer from dispersive errors in which high gradients are spread out as spurious oscillations (e.g., fig. 2.7c-d).
The numerical method used by Clawpack is an adaptively weighted combination of the second-order Lax-Wendroff method with the first-order "upwind" method (LeVeque, 2002). The upwind method is less accurate in general, but as a first-order method suffers only from dissipative errors, not dispersive errors ( fig. 2.7a-b). The Lax-Wendroff scheme is the more accurate of the two but suffers from dispersive errors in the presence of discontinuities ( fig. 2.7c-d).The second-order term in the numerical approximation controls dissipation: when this term has the correct coefficient (as in Lax-Wendroff), dissipation does not occur, but when the term has an incorrect coefficient or is neglected (as in upwind), dissipation becomes the dominant error.
Real-world pressure waveforms in gases are piecewise-continuous: they are smooth except at isolated shocks. Armed with this knowledge about the waveforms being modeled, it is apparent that preserving their shape would be best accomplished by suppressing dissipation over most of the waveform, and increasing dissipation where the gradient is very high to avoid dispersive errors. Clawpack accomplishes this by adaptive weighting of the dissipation term; the weights are provided by nonlinear functions called limiters. Limiters depend on the values along the waveform and their output is essentially a measure of wave smoothness; in smooth parts of the waveform, the limiter output suppresses the dissipative term, whereas in high-gradient areas the limiters increase the dissipative term to eliminate dispersive errors. The combination of the Lax-Wendroff and upwind methods with limiters is sometimes referred to as the Lax-Wendroff-LeVeque method, and produces much more accurate results than either constituent scheme could on its own ( fig. 2.7e-f). A wide variety of limiters have been developed, and many are included as options in Clawpack.
The CFL (Courant-Friedrichs-Lewy) number, defined in one dimension as CF L = ∆t ∆x c max (with c max being the fastest possible wave in the system) is an important quantity in any numerical model (Chung, 2002). In the Euler equations, c max is the highest value of c + |v| within the model. Model stability is partly determined by the CFL number: the condition CF L ≤ 1 is necessary (but not sufficient) for a model to be stable. Additionally, the CFL number helps determine the accuracy of numerical The first-order upwind method suffers from dissipative errors, which smooth out the signal and cause severe attenuation over time. C, D: The second-order Lax-Wendroff scheme does not suffer from major dissipative errors, but is affected by dispersive errors that cause spurious oscillations near areas with high gradients. E, F: The Lax-Wendroff-LeVeque scheme using the "MC" limiter to optimally adjust the amount of dissipation according to the local gradient. This prevents dispersive errors around shocks while minimizing dissipative errors in smooth regions of the waveform. All models shown are of the advection equation with 200 nodes (not all shown for clarity), a speed of 1 m/s, CFL number = 0.8, and periodic boundary conditions (so that the x axis represents the circumference of a circle and the waveform circles around to the starting position once per second).
schemes: in general, dissipative and dispersive errors tend to be smaller when CFL approaches but does not exceed unity.

Comparison of Finite Volume and Finite Difference Methods
By contrast, the FD method is implemented on a large number of points (termed grid nodes) distributed throughout the model region, instead of grid cells. The governing partial differential equations of the system are then discretized and implemented directly on the nodes. So, while FV works in terms of accumulated conserved quantities (stored within cells) and their fluxes (which occur from cell to cell), FD simply calculates numerical spatial derivatives at each node (thereby assuming that the spatial derivatives exist) and inserts them into the governing equations to solve for time derivatives. Consequently, the FD method does not conserve variables exactly even when the governing equations are conservation equations written as eq. 2.16; however, by not requiring conservation, it can be applied sensibly to directly model non-conserved variables like displacement or pressure.
FD methods are widely used in geophysics; for example, I use the FD method of (Kim et al., 2015) to model linear infrasound at Volcán Tungurahua in chapter 3.
However, the propensity of nonlinear systems of equations (like the Euler equations) to spontaneously form undifferentiable shocks makes FD difficult to use with such systems. Consequently, I use the FV code Clawpack to model shock waves in air.

Clawpack
Clawpack (Clawpack Development Team, 2018) is a package of codes used to solve hyperbolic systems of partial differential equations with numerical finite-volume meth-ods. Any finite-volume model must, at minimum, calculate fluxes between grid cells in each time step and update cells accordingly. Clawpack calculates fluxes between cells using Riemann solvers, which are specific to the system of equations being solved and can be implemented in different ways; many common systems of equations, including the Euler equations, are already implemented in Clawpack with one or more Riemann solvers. The work presented in this chapter uses the PyClaw branch of Clawpack; PyClaw is a Python-based implementation of the methods of Clawpack and includes a convenient user interface to configure a model run, manage variables in the model, and output data (Ketcheson et al., 2012). I also use the parallelization methods in the PetClaw (Alghamdi et al., 2011) branch of Clawpack, which enable large models to be run quickly on a high-performance computing cluster.

Assumptions of Model Simulations
I used Clawpack to simulate pressure disturbances following the Euler equations (eqs. 1.6-1.8) with equation of state given by eq. 1.9. In doing so, I assume the medium to consist of nonconductive, inviscid ideal gas, with negligible effects from gravity at these small length and time scales. Although the Euler equations do not explicitly contain temperature, I refer to gas temperature in this chapter as an intuitive proxy for density and sound speed.
The only relevant properties of the gas that depend on its chemical composition are its specific heat ratio γ = Cp C V and its gas constant R = p ρT (although, like temperature, R does not explicitly appear in the Euler equations). I assume that the propagation medium is solely atmospheric air (γ = 1.4, R = 287 J kg • C ). I assume the gas to be homogeneous in composition but variable in initial pressure and temperature. The initial condition therefore consists of a small region of high-temperature, high-pressure gas otherwise contained in air under ordinary sea-level atmospheric conditions (analogous to compressed volcanic gas suddenly exposed to the atmosphere in an instantaneous opening of the vent). Because the initial condition represents the state of the vent at the precise moment of opening, all gas is initially stationary.
Explosive yield W is a key parameter in these models. Yield can be quantified by calculating the total initial energy over the compressed gas region and subtracting the final energy in that region. This problem can be simplified with the constraints that all explosions tested in this project have zero initial fluid velocity, and that the model volume must eventually reach equilibrium at normal atmospheric pressure with zero velocity. From eq. 1.9, volumetric energy density in stationary gas is E = p/(γ−1), so total energy of a spherical explosion of radius r is W = 4 3 πr 3 pexp−1Atm γ−1 . Two explosions with equal energy but different pressures must therefore have different source radii as well. According to the ideal gas law (which can be written as pV = mRT ), total energy could equivalently be formulated as being proportional to the product of temperature T and mass m instead of pressure p and volume V .

Model Run Configuration
With PyClaw, model run settings are set in a Python script (initially provided with PyClaw as examples, but heavily modified by the user for each simulation). The basic parameters to be set are properties of the grid (number of spatial dimensions, number of nodes in each dimension, time and space intervals, and desired CFL number), the governing equations of the model and their numerical implementation (i.e., the Riemann solver, source terms, model order, and limiter), initial condition, and boundary conditions. Because the configuration script is written in Python, it is straightfor-ward to write flexible configurations that can be combined with other scripts to launch many different model runs simultaneously-a convenient ability when working with a high-performance computing cluster.

Model Results and Discussion
I ran two sets of explosion models. The first set included free-air explosions of spheres of hot compressed gas in order to determine the effects of initial pressure, temperature, and radius on recorded signals. The second set tested explosions of hot compressed gas within the conduit of an idealized volcanic edifice in order to determine effects of initial pressure and conduit geometry on recordings. The simulations discussed in this chapter are summarized in tables 2.4 and 2.5.

Free-Air Explosions: Waveform Parameters and Scaling Laws
I first ran a set of Clawpack models of explosions in the free atmosphere (table 2.4).
Each model simulated an initial condition including a sphere of high-pressure, hightemperature air surrounded by ambient atmospheric air (101.325 kP a, 15 • C for consistency with the reference explosion in (Kinney & Graham, 1985)). Models were run on a cylindrical coordinate system assuming axial symmetry ( ∂ ∂θ = 0), so the model grid included two spatial dimensions (z and r). Although this scenario is spherically symmetric as well and could therefore be simulated more efficiently in just one spatial dimension (r), cylindrical coordinates were used for consistency with later models to simulate axisymmetric idealized volcanoes.
The free-air explosion models tested and extended a basic principle of explosion study: the scalability of explosion waveforms according to explosive energy. The first test was modeling two explosions differing only in the initial compressed gas radius (with identical source pressures and temperatures). Resulting waveforms were found Subsequent tests explored the effects of source pressure and temperature on pressure waveforms in explosions of equal energy. The scaling laws eq. 2.8-2.9, which can be shown to be exact for energy differences due to changes in source radius, are not exact when energy changes due to a change in source pressure. To test this, I modeled several explosions with different pressures but equal energy and temperature. The resulting overpressure decay curves show that wave amplitudes can differ strongly vary near the source; however, the waves quickly converge to a common decay curve while overpressure is still high ( fig. 2.9). Similar results are found for the positive phase duration: all waves differ near the source, but quickly converge to equal durations ( fig. 2.10). Based solely on the distance, overpressure, and positive phase duration recorded far from a free-air explosion, it would be easy to determine explosive energy, but not explosive pressure.   Pressure waveforms from equal-energy explosions can differ visibly in shape in the trailing part of the waveform. Waveform differences for equal-energy, equaltemperature, differing-pressure explosions can be seen but are rather modest ( fig.   2.11). Temperature differences produce somewhat more visible waveform differences in the strength and timing of the trailing shock: lower temperatures result in stronger trailing shocks that arrive later ( fig. 2.12). However, in realistic field conditions in which scattering by irregular topography is significant, such small effects on the trailing shock might be difficult to obseprve due to multipathing of the initial shock.
The effect of source temperature on the properties of the trailing shock is probably a consequence of source temperature controlling sound speed of the decompressing gas. The rarefaction wave that propagates into the compressed gas is limited by the sound speed: when the sound speed is low, it takes longer for the rarefaction wave to reach the center of the gas, meaning that the resulting over-expansion, implosion, and secondary explosion will take longer, resulting in a delayed secondary shock. By contrast, an explosive gas with very high temperature and sound speed (e.g., the explosive products of TNT decomposition) can reach equilibrium so rapidly so that the trailing shock can be indistinguishable from the initial shock.
This model is consistent with the observation in fig. 2.11 that the trailing shock timing is later for the higher-pressure explosion than for the lower-pressure explo- ) to the initially compressed air, it is clear that a stronger drop in pressure corresponds to a lower post-shock temperature (and therefore lower post-shock sound speed). Therefore, for explosions of equal energy and initial temperature, higherpressure explosions will have later-arriving secondary shocks than lower-pressure ex-

plosions.
Additionally, as shown in section 2.4.2, secondary shocks are not observed when waves must propagate around a crater rim. This is noteworthy, as trailing shocks are typically not observed at volcanoes, probably due to this topographic effect.

Effects of vent geometry
The second set of simulations tested different scenarios of initial gas pressure and vent geometry. The volcanic topography was idealized as consisting of a conical edifice (with constant slope of 30 • ) containing a crater (with radius 100 m and constant slope of 45 • ) and a vertical-walled conduit with radius 50 m ( fig. 2.13b). The conduit consisted of a cylinder of variable height initially filled with hot compressed gas. The  Unlike the isotropic free-air explosions discussed earlier, the vent geometry of the models with topography (and of most true volcanoes) focuses waves upward.
Therefore, in models with topography, recordings must depend on elevation angle from the vent in addition to distance and time. In the following figures, I therefore plot pressure fields and infrasound recordings that would be expected directly above the vent and along the surface of the volcano. significantly: the amplitude of the surface pressure profile is lower and its wavelength is longer than the vertical pressure profile (both which can be expected for a wave that diffracts around a corner). Differences in the structure of the rarefaction are also apparent.
The case of fixed vent geometry with varying source pressures is explored in fig.   2.14. Source pressure has a clear effect on waveforms: explosions with higher source pressures produce waves that arrive earlier with higher amplitudes and longer positive phase durations. However, when recordings are scaled in time and taken from constant scaled distances instead of fixed points, these differences are reduced. In particular, scaling almost eliminates differences among recordings taken above the vent (fig.

2.14d).
Conversely, fig. 2.15 addresses the case of fixed source pressure (0.25 MPa) and variable initial gas height. Again, differences among these signals are apparent: larger explosions produce signals with higher amplitudes and longer wavelengths (just like free-air explosions (e.g., fig. 2.8). Important differences can be seen between the raw signals recorded on the surface and above the vent. The effect of supersonic propagation on arrival times is clear in the signals recorded above the vent but not on the surface, perhaps because of the generally higher amplitudes for the upwardpropagating wave. Additionally, although all signals recorded above the vent have  .14: Simulations of explosions with identical geometry (initial gas height of 20 m) and different source pressures. A: waveforms recorded along the surface with radial distance of 1 km. B: waveforms recorded 1.3 km directly above the vent. C: scaled waveforms recorded along the surface at a dimensionless distance of 5. D: scaled waveforms recorded at a dimensionless distance of 6 above the vent. Explosions with higher source pressures produce waves with higher amplitudes, longer positive phases, and earlier arrivals; however, these effects diminish when scaling laws are applied, particularly for recordings above the vent.
higher amplitudes than their corresponding signals recorded on the surface, the ratio between the amplitude varies and the smallest difference between the amplitudes is seen in the largest explosion, possibly due to their longer wavelength being attenuated less by diffraction. Finally, similar to the fixed geometry, variable pressure case ( fig.   2.14), application of a scaling relationship eliminates much of the differences between the signals from different explosions, especially for the signals recorded above the vent.

Discussion of Modeling Results
The numerical modeling shown in this chapter provides a starting point for inferring source properties from infrasound with a nonlinear propagation history. Scaling laws, commonly applied to chemical and nuclear explosions (e.g., Kim & Rodgers, 2016), appear to be a useful means of studying infrasound from compressed gas explosions with volcano-like pressures, temperatures, and geometries. Owing to the scaling properties of the Euler equations, scaling laws are exact in the case of free-air explosions where explosive yield varies as a result of changes to source radius only: a wavefield from a small explosion can be scaled to match a wavefield from a larger explosion exactly ( fig. 2.8). The scaling law is not exact but still useful when explosive pressure or temperature vary: essential waveform parameters like amplitude and positive phase duration can be matched using scaling laws, but subtle differences among waveforms remain (figs. 2.11, 2.12). Additionally, scaling relationships are useful (especially for signals recorded above the vent) but not exact when comparing explosions of different shapes or pressures within an idealized volcano topography (figs. 2.14, 2.15). Diffraction around the crater rim significantly affects surface-level infrasound recordings, suggesting that it may be advantageous to record infrasound MPa) and different gas thicknesses. A: waveforms recorded along the surface with radial distance of 1 km. B: waveforms recorded 1.3 km directly above the vent. C: scaled waveforms recorded along the surface at a dimensionless distance of 5. D: scaled waveforms recorded at a dimensionless distance of 6 above the vent. Explosions with greater explosive gas volumes produce waves with higher amplitude, longer wavelengths, and less anisotropy in amplitude. Again, scaling reduces this differences, especially for waves recorded directly above the vent. Waves from larger explosions also arrive earlier when recorded above the vent but not when recorded along the surface, perhaps because upgoing waves are not attenuated by diffraction and therefore propagate nonlinearly longer.
in field sites where diffraction is a less significant process (e.g., in the air above the vent or on the crater rim).
These models show that infrasound interpreted with scaling relationships can be a viable means of estimating explosive energy from compressed gas explosions at typical volcanic pressures. However, explosive energy itself is the product of quantities that may be more interesting to volcanologists (W = (pexp−patm)V Without additional constraints, it would be difficult to extend infrasound waveform interpretation to allow estimation of initial pressure, volume, temperature, or mass: changes in such properties result in only minor changes to waveforms in equal-energy explosions, which would be hard to interpret in the presence of minor uncorrected topographic scattering.
These models have limitations that keep them from being perfectly analogous to volcanoes. First, all of these models simulated explosions of hot compressed air into normal atmospheric air, whereas volcanic explosions are driven by variable mixtures of gases dominated by H 2 O, CO 2 , and SO 2 along with tephra of a range of sizes.
Even without tephra, these three gas species differ from air in their specific heat ratio (γ ≈ 1.3, vs. γ air ≈ 1.4), meaning that the thermodynamics of such gases differs from air. The behavior of H 2 O is further complicated by the possibility of phase changes.
Including very fine tephra (so that tephra particles are in mechanical and thermal equilibrium with the gas) enables the mix to be treated as a pseudogas; coarser tephra must be treated as a separate phase (Pelanti & LeVeque, 2006). Besides having complicated physics, each of these factors adds new free variables to be considered. Further, the model including volcanic topography cannot test source pressures higher than about 1 MPa due to a common tendency of models of diffracting shock waves to spuriously calculate negative densities for stronger shocks, causing the model to crash (a problem known as the "backward-facing step"). All of these are significant computational challenges that will require extension of existing CLAWPACK code to accommodate.

Conclusion
Volcanic explosions can produce waves that initially propagate nonlinearly, including supersonic propagation, changes in waveform shape, and decay. Nonlinear processes cannot be represented accurately by common analysis techniques that include convolutions, and nonlinear pressure waves must be modeled using the compressible Euler equations (which are difficult to model) instead of the acoustic wave equation (which is easy to model). Several projects have attempted to clarify volcanic shock propagation, but we still lack an explicit means of relating infrasound recordings to source physics when the wave's propagation history is nonlinear.
However, various tools for studying volcanic shocks can be adapted from past work on shock waves and anthropogenic explosions. The Rankine-Hugoniot relations, for example, establish the relationship between pressure ratio across a shock, the speed of the shock, and other fluid-dynamic properties of interest. These relations can be used in a volcanic context to estimate pressure from wave speed (if, for example, the wave speed is tightly constrained at a supersonic speed by video) or vice versa (if the pressure ratio is already known from amplitude recordings and the propagation speed is needed). A very energetic eruption is needed to create noticeably supersonic waves far from the source.
Additionally, scaling relations have an important role to play in inferences of volcanic explosions. Scaling relations show that the primary control of waveform characteristics is the distance from the source divided by a scaling parameter dependent on explosive energy. In this chapter, scaling relations were tested in the context of free-air explosions and explosions with realistic volcanic topography, and found to be useful in matching signals from different explosions (and exact in matching signals whose sources differ only in size). recorded exceptionally high-amplitude, long-period infrasound (1600 Pascals peakto-peak amplitude, 5.5-second period) on sensors within 2 km of the vent alongside electromagnetic signals from volcanic lightning serendipitously captured as interference. This explosion was one of Tungurahua's most powerful vulcanian eruptions since recent activity began in 1999, and its acoustic wave is among the most powerful volcanic infrasound ever recorded anywhere. We use these data to quantify erupted volume from the main explosion and to classify post-explosive degassing into distinct emission styles. Additionally, we demonstrate a highly effective method of recording lightning-related electromagnetic signals alongside infrasound. Detailed chronologies of powerful vulcanian eruptions are rare; this study demonstrates that diverse eruptive processes can occur in such eruptions and that near-vent infrasound and electromagnetic data can elucidate them.

Introduction
Infrasound is a valuable tool for all-weather monitoring of erupting volcanoes.
Infrasound source processes can often be resolved in detail because path effects on infrasound near the vent are often small (Fee & Garces, 2007;Johnson & Lees, 2010) or predictable and straightforward to correct (Kim et al., 2015). Local infrasound is therefore useful in tracking the style   Several volcanic processes produce distinct types of infrasound. Some of the most energetic waves correspond to discrete vulcanian (Iguchi et al., 2008) and strombolian (Gerst et al., 2013) explosions, in which explosive gas release produces brief infrasound pulses a few seconds in duration that consist of a compression, rarefaction, and coda.
The most powerful explosions generate nonlinear shock waves (Morrissey & Chouet, 1997) that decay into linear infrasound.
Other volcanic sources can generate continuous, long-duration infrasound tremor.

Field Data Collection
Our installation included two stations 1860 m and 3160 m north of the vent (Fig. 3.1), each with three infrasonic microphones with flat responses above 0.01 Hz . Lab tests on these microphone types showed negligible nonlinearity at the excess pressures measured in this study. Infrasound was logged with a RefTek RT-130 datalogger at station HIGH and a DataCube-3 at station LOW; the sample rate at these stations was 1000 Hz and 100 Hz respectively. Both data loggers recorded at 24-bit resolution with GPS timing. Sensor cables were shielded at station LOW but not at station HIGH, probably contributing to lightning-related electromagnetic interference in recordings from station HIGH (discussed in sections 3.3, 4.4, and text S2).
Terrain and equipment constraints prevented the installation of triangular arrays, so linear arrays were deployed with 15-m spacing between sensors. Although lin-ear arrays cannot be used to calculate a unique backazimuth or incidence angle of incoming waves, the spatial separation of the array elements permits us to calculate coherence among sensors, a variable we consider when analyzing tremor (section 3.4).

Numerical infrasound modeling
We The intrinsic sound speed structure of the atmosphere was calculated from the temperature profile, and the density structure was calculated using the temperature and pressure profiles ( fig. 3.3). The finite-difference method we used does not accept wind as a parameter; however, because our arrays were installed at approximately the same azimuth from the vent (∼ 17.5 • ), we used wind profiles along with the intrinsic sound speed to calculate effective sound speed along that azimuth, and used the effective sound speed as the numerical model's sound speed profile.
Timing and relative amplitudes of finite-difference results are consistent with the stations' distance to the source. However, the modeled pressure traces are compli-cated somewhat; although the signals lack prominent secondary arrivals, they are considerably longer and more complex than the 5-Hz Blackman-Harris source time function. At both stations, the main arrival lasts approximately 1 s and is followed by a coda lasting another 2 s. This complexity presumably originates as echoes and diffraction in the crater and as waves propagate over the rough volcanic topography ( fig. 3.4).

Median Filter
Infrasound recordings from station HIGH contain lightning-related glitches that need to be removed. We used a median filter for the task ( fig. 3.2). Median filters approximate linear running average filters over non-glitchy data (in which the median is approximately equal to the mean) but are highly effective at removing glitches (in which the median is very different from the mean).
Median filters have just one parameter to adjust: the window length. Longer window lengths result in a lower corner frequency of the filter and more effective glitch suppression; shorter windows have a higher corner frequency but are not as effective at removing glitches. Fig. 3.2 shows two median filters with different window lengths.
Both are effective at removing glitches. We selected the longer window length (0.061 s, 61 samples) which approximates a linear low-pass filter with corner frequency (-3 dB) of 7.3 Hz. We selected this filter because the higher frequencies that are affected by the longer window are not important to our analyses. In general, the negative effects of median filters can be mitigated or avoided by recording at a high sample rate.

Analysis
Infrasound analysis enables us to quantify eruptive activity by the amount of air it displaces (summarized in table S1). We note that some of these analyses address details of waveforms that can be altered by path effects; therefore, recording infrasound near the vent where path effects are weak is essential.
Infrasonic pressure can be used to calculate erupted volume during discrete events.
For linear acoustic waves emanating from an isotropic source, pressure recorded at the microphone corresponds to displaced air at the vent as where q(t) is mass flow rate of displaced air, H(t) is the inverse of the Green's function from the source to the receiver, *' denotes convolution, and p(t) is the infrasound time series in pressure units. When a signal is recorded at multiple sites with different Green's functions, an overdetermined linear system can be constructed and solved: where H ij is the ordinary-least-squares generalized inverse of the Green's functions, and p j is the recorded data from all receivers (using the Einstein summation convention).
In the particular case when path effects (such as from topography, atmospheric heteogeneity, and attenuation) are negligible, equation (3.1) can be simplified as where r is the distance between vent and microphone and Ω is the solid angle (about 9.08 steradians at Tungurahua) subtended by the atmosphere around the vent (Lighthill, 1978). Otherwise, equation (3.2) must be used instead with Green's functions calculated using numeric models.
Eqs. 1-3 make four assumptions (Johnson et al., 2003). First, any instrument response, path or site effects, or radiation pattern anisotropy must be accounted for in H (and must be negligible if using equation 3.3). We calculate Green's functions using finite-difference modeling with appropriate topography and atmospheric structure (figures 3.3, 3.4).
Second, long-period noise (such as from wind) and instrument drift (small trends in recordings not representative of actual pressure changes) must be negligible. Such noise is common in infrasound recordings and is problematic because trends and low frequencies are magnified by inversion, sometimes causing volume flow to be non-zero long after the end of the signal. Johnson & Miller (2014) solved this problem by estimating flow rate using equation (3.3) and detrending to make infrasound pressure and flow rate zero at the end of the signal. We use a slightly different detrending method because we use Green's functions from numeric models instead of equation 3.3. Using linear inversion, we find an optimal drift that minimizes the L2 norm of estimated flow rate:  where m is the slope of the instrument drift and t j is time (using the summation convention). Apart from instrument drift, signal-to-noise ratio is high during the explosion (Fig. 3.5).
Third, the source dimensions must be small compared to acoustic wavelengths to justify a point source approximation. With vent dimensions of tens of meters compared to dominant wavelengths exceeding 1 km, this assumption is valid.
Fourth, pressure perturbations must be small relative to ambient pressure to justify the use of linear acoustic theory. This is easily satisfied for the precursory vent uplift but possibly not for the explosion. The wave from the explosion exceeded 1200 Pa peak pressure at 2158 m, so the wave's excess pressure probably reached 10% of ambient atmospheric pressure within 400 m of the vent. Nonlinear effects in this region could have caused waveform change and decay, which equations 3.1-3.4 cannot account for. However, significantly nonlinear propagation typically forms an abrupt signal onset, which we do not observe in this signal. Although we suggest that nonlinear effects on this wave were probably weak, we must consider our volume calculation to be conservative.

Main Vulcanian Blast
The eruption began abruptly with a major explosion that produced the high-amplitude infrasound signal (Fig. 3.5). Like waves from many explosions, the waveform is dominated by a strong compression followed by a longer, weaker rarefaction (Morrissey & Chouet, 1997). However, unlike many explosion waveforms, the rise from ambient to peak pressure occurs in several distinct steps; pressure rises quickly at the onset of the waveform and 0.6 seconds, 1 second, and 1.4 seconds later (Fig. 3.5).
Two of the three channels at station LOW clipped during the peak of the main blast arrival; consequently, data from those two microphones are excluded from the blast analysis. Only the main blast arrival was clipped in these two channels, so we did consider them in other analyses. The remaining channel at station LOW did not clip due to its lower sensitivity, and the high input voltage range of the RT-130 logger prevented clipping on any channel at station HIGH.
We apply equation (3.4) to quantify erupted volume using data from the four microphones that did not clip. We first estimate the volumetric eruption rate sensorby-sensor and find that one sensor (HIGH-3) disagrees with the others (Figs. 3.5e-f).
This outlier is attributed to an uncorrected instrument drift or other long-period noise, and we omit it from subsequent calculations. The remaining three sensors agree closely with each other (±10% of cumulative flow). We invert these three records jointly for erupted volume, finding peak flow rate of 2.39 × 10 7 m 3 /s and cumulative volume of 5.31 × 10 8 m 3 .

Precursory Uplift
The first infrasound signal from this eruption was a relatively small emergent pressure increase (beginning at 11:46:38 UTC at the nearest station) that lasted about 0.71 s before being obscured by the main blast arrival (Fig. 3.5). No coherent infrasound was recorded before the onset of the emergent pressure rise.
Inverting the brief, low-amplitude precursor using Green's functions from numeric models would be problematic because of possible acausal contamination from the blast wave. Instead, we invert the precursor using equation 3.3 (which is strictly causal) and correct for path effects by weighting each trace by the amplitude of the corresponding Green's function. Drift removal is performed by subtracting the trends found in the 14 seconds before the precursor began. Signal-to-noise ratio is not high during the precursor, so weighted traces are stacked to reduce noise.
We estimate the volume of displaced atmosphere during this uplift as 2.08×10 4 m 3 using equation 3.3 (Fig. 3.5). Following interpretations of similar precursory pressure increases at Santiaguito (Johnson & Lees, 2010), Sakurajima (Yokoo et al., 2009), and Suwanosejima (Yokoo & Iguchi, 2010), we attribute this signal to vent surface uplift before gas escapes. The area of the vent that deformed during this phase is unknown, so we cannot calculate the height of the uplift.

Electrical Activity
Infrasound from station HIGH contains glitches (brief spikes in the time series) throughout the first 20 minutes of the eruption. Glitches do not appear at station LOW, possibly because of differences between the two loggers' locations, antialiasing filters, or cables (which were shielded at LOW but not at HIGH). These glitches include a one-sample voltage spike preceded and followed by brief oscillations (presumably induced by the data logger's antialiasing filter); they appear on all channels within the same sample interval (0.001 seconds) with distinct waveforms ( fig. 3.2).
Such characteristics are unlikely for acoustic waves, but are typical of glitches observed to coincide with lightning strikes during thunderstorms, which are interpreted as electromagnetic interference from radio waves generated during lightning strikes (Anderson et al., 2014). Volcanic lightning is common in explosive eruptions and forms as a result of charging due to violent interactions between ash particles or, if present, hydrometeors (Behnke et al., 2013;Cimarelli et al., 2014).
Glitches are considered noise for the infrasound analysis and are therefore removed by a median filter before waveform inversion (text S2). Although glitch elimination provides cleaned-up acoustic signals to analyze, we also analyze the timing of the glitches because they indicate lightning activity during the eruption. Lightning is absent from 0-25 s, frequent from 25-250 s, sparse from 250-1250 s, and absent after 1250 s (Fig. 3.6). In total, we identify 119 lightning events, of which 93 occurred within 300 s of the explosion. For comparison, the World Wide Lightning Location Network (WWLLN, a network that detects lightning globally with low detection efficiency) detected three events, all of which appear in our data.

Tremor
The main vulcanian explosion is followed by a period of infrasonic tremor. We subdivide the tremor into types 1a, 1b, and 2 by waveform shape, amplitude, and semblance ( Fig. 3.6). The first tremor period (60-1200 s) alternates between tremor 1a and 1b, and the second tremor period (1200-3000 s) consists of type 2 only. The following analyses are summarized in table S2. All analyses are done with data from channel HIGH-2.
Waveform shape varies systematically during the long period of tremor following the explosion, with some periods composed of strongly asymmetric pulses (higheramplitude compressions than rarefactions) and other periods fluctuating more evenly about zero. To quantify pulse asymmetry over long time scales, we plot Pearson's moment coefficient for skewness for moving 50-s windows. Skewness is defined as where µ 3 is skewness, µ is the mean, σ is standard deviation, and N is the number of samples in a window. This measure quantifies the degree to which tremor is dominated by strongly asymmetric pulses . Additionally, we calculate semblance as an indicator of wave coherence. For a window of n samples, the semblance S is defined as where p j is the infrasound recorded at microphone j, t i is the time step i, and δt j is the time shift at microphone j for the apparent velocity. For each window, we test a range of apparent velocities and use the maximum semblance value obtained.
Tremor 1a and 1b both feature continuous, stationary waveforms. However, they differ from each other in skewness, amplitude, and semblance. Skewness of tremor 1b fluctuates around zero, while tremor 1a ranges from about zero to one. Further, both semblance and amplitude are greater in tremor 1a than tremor 1b.
Tremor 2 differs from tremors 1a and 1b mainly by waveform shape. Unlike the continuous waveforms of tremor 1a and 1b, tremor 2 waveforms consist of closely spaced asymmetric pulses separated by irregular time intervals. Each pulse includes a strong, short-duration compression followed by a weak, longer-duration rarefaction; amplitudes of compressions are 2-3 times those of rarefactions. Correspondingly, skewness rises sharply at the onset of type 2 tremor and remains high throughout that period, demonstrating that the visible differences in waveform shape constitute a large-scale structural difference between these tremor types. Zero-to-peak amplitudes of tremor 2 are similar to those of tremor 1a, but peak-to-peak and root-mean-squared amplitudes are lower in tremor 2; this discrepancy results from the high skewness of tremor 2.

Discussion
We divide the volcano's activity into three periods, including the main vulcanian explosion and two distinct degassing phases.

Main Vulcanian Explosion
Before the explosion, the vent was sealed and contained a large quantity of pressurized gas. Surface activity began with a rapid uplift of the vent, creating a precursory acoustic pulse lasting at least 0.71 s and displacing at least 20800 m3 of air. Because the shock wave produced by gas release may have propagated supersonically and partially overrun this precursor, the uplift might be underestimated.
Explosive gas and tephra ejection followed the vent uplift. The stepped rise of the wave onset in the infrasound data ( Fig. 3.5) suggests that gas release occurred in pulses, probably due to incremental opening of the conduit or successive tapping of deeper gas-charged sections of the magma column. Altogether, the total erupted volume from the main explosion is estimated as 5.31×10 8 m 3 including gas and tephra . If, contrary to our assumptions, nonlinear propagation effects were significant, these values would be underestimates.
To provide context for the scale of the 14 July 2013 eruption, we compare it to the two largest explosions from a recent period of explosive activity at Sakurajima volcano, Japan (explosions 3 and 5 of Johnson & Miller (2014)). Infrasound analysis indicated that the 2013 Sakurajima explosions erupted 8.3 and 8.4 ×10 6 m 3 of volcanic gas and tephra, a factor of 63 less than the volume erupted at Tungurahua. Compared to these eruptions, Tungurahua's blast wave had a longer compression duration (2.4 s compared to about 1 s at Sakurajima) and a higher peak-to-peak reduced amplitude (3.4 ×10 6 P a·m compared to 3.9×10 5 and 7.4×10 5 P a·m at Sakurajima). Therefore, the Tungurahua explosion was considerably larger than these Sakurajima explosions by several measures, especially by infrasound-inferred erupted volume. By infrasound amplitude, this explosion was also more powerful than any other discrete explosion at Tungurahua since the current monitoring network was installed in 2006. In particular, the peak-to-peak amplitude of the July 2013 explosion exceeds amplitudes of later vulcanian eruptions in October 2013 and February 2014 by factors of 1.7 and 2.9, respectively (Hall et al., 2015).
In this explosion, the vent opened in a complex 2-second process including preexplosive vent deformation and a series of emission pulses. Similar multi-pulse superposition of blasts is evident, though less pronounced, in explosions at Sakurajima Johnson & Miller, 2014), and similar pre-explosive deformation has been observed at volcanoes including Sakurajima (Yokoo et al., 2009), and Suwanosejima (Yokoo & Iguchi, 2010). These observations demonstrate the complexity of vent processes that can initiate powerful explosions: vent opening can be a series of events rather than one single failure. By comparison, the explosive mechanism was much simpler (single explosive pulse, no pre-explosive deformation) in Tungurahua's small explosions in the weeks following 14 July 2013.

Continuous Tephra-rich Degassing
Infrasonic tremor oscillating between types 1a and 1b follows the explosion and is accompanied by sporadic electrical discharge. The lack of discrete waveforms within the tremor indicate a relatively continuous emission process, although the variation in amplitude between these tremor types shows that emission vigor fluctuated over time scales of tens to hundreds of seconds.
During tremor 1a periods, the high semblance among sensors corresponds to high signal-to-noise ratio and a single dominant acoustic source. As a result, we consider this signal to indicate continuous magma fragmentation at the vent, either by downward propagation of a fragmentation wave (Alidibirov, 1994) or continuous ascent and bursting of bubbles (Ulivieri et al., 2013).
We interpret the low-semblance signals of tremor 1b to be dominated by incoherent pressure variations that mask eruptive signals during periods of reduced emission vigor. This incoherence probably arises from severe wind noise, possibly driven by updrafts near the vent, as station HIGH was installed in a sparsely vegetated site with no wind reduction system. Other possible sources of incoherent noise may include PDCs that passed within hundreds of meters of station HIGH (Hall et al., 2015) and falling ejecta.
Electrical discharge began before the first tremor period and continued throughout. The incidence of lightning-linked glitches agrees with the pattern seen in recent studies at Augustine (Thomas et al., 2007), Redoubt (Behnke et al., 2013), and Eyjafjallajokull volcanoes (Behnke & McNutt, 2014). These studies describe explosive eruptions in which electrical discharge started near the vent soon after explosive activity began. Two styles of discharge occur during this period: vent discharge (small events occurring continuously at the vent), and near-vent lightning (higher-power discrete lightning flashes near the vent). A third type of discharge, referred to as plume lightning, occurs when a plume is present; plume lightning occurs in long, powerful, discrete flashes similar to thunderstorm lightning (Behnke et al., 2013).
Similarly, we observe a period between about 25-250 s after the eruption onset with frequent discharges (around several tens of discharges per minute). The beginning of this phase is abrupt and no discharges are observed earlier. These discharges are probably near-vent lightning; our instrumentation is not designed for sensitivity to electrical activity and is therefore unlikely to record the lower-power vent discharge.
After about 250 s, electrical discharge becomes much less frequent (around 1-2 discharges per minute). Because we cannot locate these events, we cannot classify them between plume lightning and near-vent lightning. The period of sparse lightning ends around 1250 s, and no further discharges occur after that time.

Pulsed Degassing
The properties of the tremor signal change abruptly around 1200 s. This new signal, labeled tremor 2, is dominated by discrete pulses separated by irregular time intervals.
Pulsed degassing continues for 1800 seconds before diminishing.
Electrical discharge also ceases around the beginning of the pulsed degassing phase. Two interpretations of plume activity could explain the absence of discharges.
The first possibility is that plume electrification was driven by volcanic emissions during the first tremor period, so lightning diminished when the style of emission changed. The second possibility is that plume electrification resulted from an initial intense venting of gas and ashnot by emissions during the first tremor periodand that by coincidence the plume electrical activity stopped around the transition to tremor 2. In either case, emission during the tremor 2 period was insufficient to create a column or cloud with electrical activity.
We consider two potential sources for tremor 2, one being small, repeated bursts at the vent. This mechanism is suggested by the resemblance between these pulses and explosion waveforms (sharp, brief compressions followed by longer, weaker rarefactions). These events could be bubbles ascending to the magma surface and bursting (e.g., Ulivieri et al., 2013), or cycles in which permeable crack networks in the magma pressurize, burst, and re-seal.
These waveforms also resemble crackle, a type of signal produced by supersonic jetting. Waveforms similar to tremor 2 were recorded at regional distances during the June 2011 eruption of Nabro volcano (Eritrea) and Stromboli volcano (Italy), and were attributed to supersonic gas jetting because of their similarity to noise from rocket and jet engines Goto et al., 2014). Crackle was associated with the remarkably ash-poor, gas-rich character of the Nabro and Stromboli emissions.
We speculate that ash-poor jetting could potentially explain both the absence of electrical activity and this tremor.

Recording Volcanic Lightning Alongside Infrasound
Our method of lightning glitch loggingconnecting infrasound microphones to a data logger with long unshielded cableswas discovered by accident in this campaign but will be useful in future infrasound projects. If a local infrasound station is already being installed, this method yields lightning event times with no additional equipment, expense, or installation time, and minimal extra data processing, and with higher detection efficiency than is possible with global networks like WWLLN. The timing and event frequency of volcanic lightningwhich is driven by conditions in the plumecan therefore serve as an easily obtained but valuable complement to existing monitoring data (which typically offer little information on plume conditions in cloudy weather).
However, this method has downsides and is inappropriate in some scenarios. For example, if more detailed information (like discharge location and power) is required, the method in this paper is insufficient and a dedicated lightning monitoring system should be used (Behnke & McNutt, 2014). Additionally, because the median filter acts as a low-pass filter on ordinary data, a higher sample rate than normal may be needed to prevent the median filter from attenuating high acoustic frequencies of interest.

Conclusion
We The short duration of campaigns presents researchers with an opportunity to collect valuable data at low cost and with few workers. This is beneficial, especially to researchers with limited resources, but imposes constraints in the field: many difficult tasks must be accomplished with limited time and personnel. For example, the common need to install equipment in remote sites (for reasons like optimizing signal-to-noise ratio, or to make a network spatially comprehensive) often requires long hikes and camping, limiting the amount of equipment that can be carried as well as the time available for installation. In general, the number of stations that can be installed during a campaign is constrained by instrument portability and ease of installation, and campaigns that use many portable, easy-to-install instruments can be more effective than those that rely on few heavy, expensive systems.
Instrumentation needs during campaigns differ from those of long-term projects.
Power consumption of sensors and loggers is critical because it affects portability: for example, a single infrasound station that requires a 10-kg lead-acid battery may be less portable than several stations powered by 500-g sets of alkaline batteries.
Additionally, time constraints during campaigns mean that installation and retrieval procedures should be reduced to as few tasks as possible (ideally just site selection, note taking, and connecting or disconnecting power). Simplifying these procedures also reduces the probability of data loss occurring due to user error during potentially stressful conditions encountered in the field. For example, combining the sensor and logger into a single unit can improve the instrument's portability and ease-of-use by eliminating the need for cables (discussed in section 1.2).
Finally, cost determines the number of instruments a researcher can afford to install, but also affects the quality of data that can be recorded. Equipment in the field faces a variety of risks, including vandalism, animal attacks, and geologic hazards. (All of these have affected our instruments in past infrasound campaigns: the authors have had sensors stolen, damaged by bears, flooded, and buried under lava.) Researchers face a dilemma regarding signal quality and instrument vulnerability: the best sites are often the most vulnerable, with the area near a volcano's vent being a prime example. When weighing the benefits of better data against the risk of instrument loss, a cautious researcher would opt for riskier sites with higher data quality when instruments are cheap, and for safer sites with lower data quality when instruments are expensive. Therefore, the use of cheaper instruments enables researchers to install larger networks that include high-risk, high-reward sites.

Inadequacy of Existing Data Loggers in Infrasound Campaigns
Various high-quality seismic data loggers are often used in infrasound campaigns.
Examples include the Omnirecs Datacube (3 channels) as well as the Trimble REF TEK 130 (3-6 channels). Each of these loggers records 24 bits at a range of sample rates with multiple gain settings, which makes them compatible with a wide range of sensors beyond seismometers, including infrasound sensors. A third datalogger to consider is the Trimble REF TEK 125A, a single-channel datalogger designed to log geophones in short deployments during active-source seismic surveys. A limitation of these data loggers is their expense (at least 1000-2000 USD per channel). These costs are comparable to broadband seismometers, but exceed by an order of magnitude the cost of some common infrasound microphones (such as the infraBSU, described in ). This cost disparity limits the opportunity for large-N projects and deployments in risky or restrictive environments.
Further, multi-channel data loggers require long cables for effective array geometries, which is logistically problematic for campaign-style infrasound deployment. Cabling limits the aperture of microphone arrays and is expensive, bulky, hard to conceal, prone to pest damage, and generally inconvenient. On the other hand, arrays of single-channel sensor-loggers may be distributed with a high degree of flexibility and arbitrarily long spacing. In practice, microphone arrays using single-channel instruments are much easier to conceal and faster to install, and higher-quality scientific results can be obtained from their more flexible geometries (Christie & Campus, 2010).

Designing a Custom-built Infrasound Logger
We first decided to explore designing and building a custom infrasound logger to overcome the limitations of existing infrasound logging systems. The Gem was designed with the priorities of portability, ease of use, power draw, and low cost. The logger is integrated with a calibrated infrasound sensor ; this allowed simplifications to be made that reduced the package's cost, power draw, and size. The resulting instrument is ideal for infrasound campaigns lasting from hours to months.
The design process for the Gem has lasted around four years off-and-on at the time of this writing (Table 2) and will probably continue in the form of minor bug fixes and improving extensibility of the logger. Development began with finding and testing methods that could record powerful infrasound at 100 samples per second with high resolution. Ultimately, a 16-bit analog-to-digital converter combined with an infrasound sensor  was found to work. Subsequent functions (writing to a micro-SD card, receiving GPS data, and using an efficient power supply) were added successively, ultimately resulting in a basic infrasound logger. This first version was merely a proof of concept: signal resolution was poor due to the lack of an amplifier, assembly was tedious because a printed circuit board (PCB) was not available, and the makeshift enclosure caused loss of function due to severe overheating in the sun.
These problems were addressed in the coming months, resulting in the first numbered version (0.5). The main improvements to version 0.5 included adding a lownoise amplifier (so signals could be recorded at high resolution), designing a PCB (making the manufacturing process less tedious and error-prone), and finding a suitable enclosure that was convenient, watertight, and inexpensive. The need to learn new design skills (PCB design software for example, discussed further in section 4) and the details of unfamiliar components made progress slow. Mistakes from inexperience led to more delays: for example, omitting a protective soldermask from the PCB inadvertently made soldering much more difficult and time-consuming. Once complete, version 0.5 was shown to be capable in the field, and some units are still in use. However, it had much room for improvement: the logger had insufficient protection from power supply noise, manufacturing was difficult and slow, and care was needed during use to avoid short circuits.
Subsequent versions made small changes to hardware (interchangeable parts, reducing self-noise, and improving manufacturing methods) and firmware (improvements to user interface, adding a digital anti-aliasing filter). The most significant changes occurred in version 0.9, in which most through-hole components and breakout boards were replaced by compact, inexpensive surface-mount equivalents. Additionally, this version incorporated a custom microcontroller board with a built-in efficient power supply, as well as an analog power plane powered by a noise-reducing linear regulator. In this paper, all descriptions of the Gem's design and specifications refer to version 0.9, and firmware code and hardware designs for this version are provided in the electronic supplement to this article (codes S1-S3). Gem development was mainly a side project throughout most of the process, although synergistic application of the Gem to infrasound research projects provided essential testing opportunities (described in section 3).
We present this instrumentation to demonstrate the feasibility for individual scientists to develop their own specialized instruments. Hardware design does not necessarily require the resources and expertise of a dedicated instrumentation company.
Development of this sensor-logger package began with integration of DIY-oriented components, some of which are still used in the current version. We used free, opensource software to program the microcontroller in a common language (C++) and to design the printed circuit board. Electronic components, soldering equipment, a computer, and an inexpensive programming cable were the only supplies needed in this process. The authors are not professional engineers and we often relied on hobbyist-oriented online forums for support.

Theory of Operation
The Gem consists of an infrasound sensor , analog signal processing circuitry, microcontroller, analog-digital converter (ADC), GPS receiver, micro SD slot and card, power supply, and sensors to track battery voltage and temperature (figure 4.1). The electronics occupy two stacked circuit boards, the lower containing the microcontroller and power supply and the upper containing all other components. Elements of the popular Arduino platform (ard, n.d.) are used in this system: the microcontroller runs the Arduino bootloader, Arduino software is used to program the microcontroller, and the standard pin arrangement of the Arduino Uno is used for connecting the two boards.
The microcontroller samples infrasound via the ADC at 400 samples per second (sps), applies an anti-aliasing filter, and writes decimated output to the SD card at 100 sps. State-of-health, GPS location, and GPS time data are also collected and written to disk periodically. Data are stored as hour-long text files, with each data type (infrasound, state-of-health, and GPS) written as formatted lines. (Text is preferable to more-compact binary because text files are easy to read and troubleshoot, and disk space is not a significant constraint.) Once data acquisition begins, the Gem logs continuously until the battery dies or the user stops acquisition.
Many aspects of the Gem are user-configurable. Amplifier gain is controlled by a resistor that can be replaced by the user. Hardware components are modular, allowing straightforward changing of the pressure transducer model (which controls the dynamic range) and the capillary filter assembly (which controls the sensor's The infrasound sensor produces an analog signal that is low-pass filtered by a resistor-capacitor (RC) network, amplified by an instrumentation amplifier, and digitized by the analog-digital converter (ADC). The microcontroller collects data from the ADC, Global Positioning System (GPS), temperature sensor, and battery voltage sensor, and logs data and environmental information to the SD card. A power supply converts battery power into a steady 3.3 V that powers the digital components; to reduce noise, a linear regulator powers the analog components.
(b) Annotated photo of Gem logger electronics (without enclosure). (c) Gem logger prepared for field use in a typical watertight enclosure with rechargeable battery. Solar power is provided through a plug on the outside of the enclosure.

Signal Conditioning
Data loggers must eliminate spectral energy above the Nyquist frequency while preserving signals below the Nyquist frequency, and must digitize data with high resolution. To accomplish this, the Gem incorporates filtering stages to eliminate high frequencies and amplification to ensure high-resolution recordings. These conditioning steps are needed to avoid aliasing and to take full advantage of the dynamic range of the signal.
The sensitivity of the pressure transducer is low (20.4µV /P a), so amplification is necessary to align typical output voltages with the input range of the ADC. We use an instrumentation amplifier to apply a gain of 23.5 to the transducer output to make the combined sensitivity 478.4µV /P a. The instrumentation amplifier requires a lowimpedance reference voltage, which is provided by an operational amplifier tracking a voltage divider. Finally, the ADC samples the amplified signal at 16 bits, with a resolution of 7.8125µV (16.3 mPa) and clipping amplitude of +/-256 mV (+/-536 Pa).
The Gem includes several filters. First, clipping from long-period pressure drifts of the ambient atmosphere is prevented by a single-pole high-pass pneumatic filter based on a capillary tube and backing volume . Second, a resistor-capacitor low-pass filter suppresses high frequency (¿200 Hz) analog signals to prevent aliasing. Third, the ADC averages the input signal over the 2.1-ms sample period, which further reduces high-frequency fluctuations. These last two filters are necessary to prevent aliasing; however, their roll-off is relatively gradual, so a digital filter is also incorporated to provide a sharper transition. To do this, we oversample at 400 sps and apply a digital low-pass FIR filter with 25 sinc coefficients, eliminating frequencies above 50 Hz. This allows the signal to be downsampled and written to disk at 100 sps.
The combined effect of these filters approximates our objective of eliminating high frequencies (¿50 Hz) while preserving lower frequencies ( figure 4.2). Aliasing should be minimal except for cases with weak infrasound and extremely intense audible signals at particular frequencies. Because most volcanoes (and geophysical sources in general) emit acoustic energy mainly in the infrasound band (Johnson & Ripepe, 2011), we consider the possibility of aliasing to be remote in realistic field situations.

Real-Time Operating System
Another requirement of the Gem is to digitize data at a precise sample rate with equal intervals between samples. Although writing to the SD card usually requires little time, write times can vary and potentially exceed the 2.5 ms sample interval, which would compromise sample time accuracy. Therefore, we use a real-time operating system (Greiman, 2017) including two threads to assign tasks different priorities: a primary thread addresses tasks directly related to sampling, and a secondary thread performs all other tasks (including processing GPS strings, recording state-of-health data, and writing to disk). The primary thread interrupts the secondary thread at the beginning of each sample interval, records a sample, and then returns control to the secondary thread. This allows coordination of both high-and low-priority tasks, ensuring that sample times are correct.

GPS Timing
Precise timing is essential for infrasound data loggers and for array processing of infrasound. The Gem obtains timing information from a built-in GPS module. Every second, the GPS provides two types of data: a brief pulse marking the exact beginning of the second, and a text string including the time, date, and location. The logger writes this information to the SD card so that correct sample times can be determined during post-processing.
Each GPS, metadata, and data line includes a timer count (accurate to ±30 ppm) from the microcontroller. The data conversion software determines the relation between the timer count and GPS time and uses the relation to resample data at regular 0.01-s intervals. This procedure enables accurate sample times to be obtained without the computational expense of real-time clock drift corrections.
The GPS consumes more power than any other component in the Gem, so it is activated briefly every fifteen minutes, running long enough to record twenty time stamps. This GPS cycle is sufficient to maintain timing accuracy of ±2.5ms. Sample time accuracy was confirmed by comparison against trusted data loggers in lab and field tests (for example, fig. 4.3).

Metadata
State-of-health information is logged every 10 s. These metadata streams include battery voltage and internal temperature, along with information on memory use and processing time. Such information is useful for understanding battery life and temperature effects on batteries under realistic field conditions, as well as for debugging and investigating unusual recordings. Additionally, the user may connect up to two extra analog sensors whose output is recorded and logged as metadata. This allows users to customize the Gem to record additional environmental information appropriate for their particular projects.

Power
The power supply for the Gem must be efficient (to maximize battery life) and steady (to optimize signal quality). However, ordinary Arduino boards are not very efficient: components not needed in ordinary operation draw significant current, and power is provided by a wasteful linear regulator. To increase power efficiency, we use a custom Arduino-like microcontroller board (referred to as the Efficientino) that in-

Performance
This section evaluates the cost, ease-of-use, and accuracy and reliability of the Gem.

Cost
Costs of components depends on quantity ordered. Custom-built sensors are often produced in relatively low quantities; assuming a ten-unit build, the total component cost is less than about 200 USD per Gem.
Labor costs must also be considered during unit construction. A worker with decent soldering skills can build a Gem in 2-4 hours. At a university, student hourly labor keeps the cost low such that per-unit construction is less than 300 USD. We find that the build process helps students by developing their electronic skills and introducing them to the field of engineering and scientific measurements.
Overall, the cost of building the single-channel Gem (with built-in sensor) is less than 10% of the cost of a conventional seismic-style datalogger even without sensors.
Per-channel, the Gem is less than about 30% of the conventional logger's cost.

Field tests at volcanoes
We tested the Gem in a range of volcanic settings to determine its accuracy, reliability, and convenience in the field. Some tests were conducted in parallel with installations using commercial 24-bit dataloggers. Initial installations were carried out at two We also performed several tests in controlled environments to assess the noise level of the Gem. An appropriate testing environment should have GPS reception and low levels of noise from ambient infrasound and electrical interference, so typical lab environments are often unsuitable for this task. In these tests, multiple Gems are co-located and subjected to a common infrasound signal; the resulting recordings should agree to high precision within the frequency band of interest. The results of these tests agree with the field tests and establish a noise floor of 7.2-15.1 mPa RMS (depending on frequency band) for the most common configuration (table 1).

Volcan el Reventador
Two Gems were installed near Reventador and co-located with an infrasound microphone logged by a commercially available 24-bit Omnirecs Datacube. Although the instruments were installed by field personnel who had no previous experience with the Gem, it was reported to be easy and convenient to install and transport. Additionally, the environment at Reventador is very wet and the Gems remained watertight during the nine-month installation period.
Recovered data quality was excellent as quantified by the high degree of similarity between recordings from the Gem and Datacube (figure 4.3). In two explosions studied, correlation between Gem and Datacube recordings was 0.996 and 0.995 over the 0.5-10 Hz band. During the ten seconds preceding these two explosions, the Gem and Datacube signals differed by 0.019 Pa and 0.027 Pa RMS in the same band, again an excellent level of agreement.

Volcan Fuego
The Fuego project was conducted to test the usefulness and ease of installation of a large network of instruments during a short four-day project. A single worker carried ten Gems in addition to camping equipment, and installed the Gems alone at sites high on the volcano during one morning of work.
This test demonstrated the utility of the Gem in the field. The light weight (3.6 kg total weight of all instruments), small size, and fast installation procedure meant that little time was required to install the network. A great deal of deployment flexibility was gained by the single-channel nature of the Gem: the network geometry used at Fuego would have required hundreds of meters of cable if traditional multi-channel loggers had been used. Further, with no long cables to bury, it was easy to conceal the instruments to protect against theft and vandalism in this theft-prone site.
Signals from the network of Gems were used to invert for infrasound sources using the finite-difference models of (Kim et al., 2015) and the explosion gas flux estimate method of (Johnson & Miller, 2014). In this procedure, a finite-difference model is

Non-Volcanic Applications of the Gem
The Gem is specialized for campaign-style fieldwork on volcanoes. However, its qualities of being portable, inconspicuous, and inexpensive have made it useful for other infrasound studies, a few of which are summarized in this section.

Aerial Infrasound Recording
Infrasound is normally recorded on the ground. However, airborne infrasound has received recent attention such as the balloon-based study of (Bowman & Lees, 2015b).
In such work, payload weight and size is limited, so it is important that instruments be small and lightweight. Additionally, aerial platforms are inherently risky, making low instrument cost desirable. These demands make the Gem an attractive instrument for such applications.
Several Gem loggers have now flown on high-altitude balloons, recording infrasound tens of km above the surface of the Earth. Two types of balloons have been tested as platforms for aerial infrasound: solar hot air balloons, and helium balloons.
Recently, a solar balloon was designed that can be built and launched anywhere by a single worker, requiring several hours of preparation and costing a few tens of USD in materials (Bowman et al., 2015). Although the balloon can reach altitudes of 16-22 km and can fly as long as sunlight is available, it has an extremely low lift-to-size ratio and therefore cannot carry heavy payloads. In fact, the Gem is the only known infrasound recording system that it can carry safely.
In a recent project, a solar balloon-borne Gem at an altitude of 16 km recorded a controlled explosion on the ground at a range of 320 km ( fig. 4.5), demonstrating the utility of balloon-borne Gems for recording infrasound in order to locate sources or probe atmospheric structure. The weight of the Gem was reduced to 77 g by removing it from the standard enclosure, and a lightweight 9-V lithium battery was used to power it. The amplifier gain was increased by changing the resistance of the gain-setting resistor, resulting in resolution of 4.87 mPa/count; this helped improve signal fidelity when recording weak infrasound in the stratosphere. Background noise levels on the Gem logger were comparable to those recorded on the International Monitoring System infrasound network ( fig. 4.6). The 0.2 Hz ocean microbarom spectral peak was present during the flight, but ground stations in the balloon's flight path did not record this signal. The lack of a peak on the ground is consistent with high wind noise levels during the day (Roger et al., 2005), an issue that balloon-borne sensors do not seem to have (Bowman et al., 2017).
Gem loggers are also useful as piggyback payloads that can take advantage of existing balloon campaigns to access regions lacking in traditional infrasound coverage, like the open ocean and the poles. These large helium balloons are typically carry instruments such as telescopes that weigh thousands of kilograms; a Gem logger's weight is trivial by comparison. However, the scientific impact of recording infrasound in the stratosphere is significant (e.g., Bowman et al., 2017), so balloon science teams are often willing to include such lightweight additions to their main payload.
Gems have been tested in this capacity already: they were present on the 2015 High Altitude Student Platform flight (Guzik et al., 2008;Bowman & Lees, 2016), where they successfully operated during a day/night cycle at an altitude of 37 km.   (Brown et al., 2014). Gem power spectra as recorded (unscaled) and adjusted for decreased acoustic impedance at altitude (scaled) are presented (Bowman et al., 2017).

River Infrasound
The Gem was recently used to record infrasound from hydraulic jumps on rivers. We studied several features in the Boise and Payette rivers, Idaho, USA. As an example, we discuss a campaign at Big Falls, a stepped series of hydraulic jumps on the South Fork of the Payette river (Ronan et al., 2017).
The purpose of this project was to record river infrasound through a flood cycle to track hydrodynamic changes in the rapid and their relationship to river discharge.
This was accomplished by installing an instrument network around the rapid. This site was difficult to access-in particular, the left bank could only be reached by whitewater kayak-and had a high risk of vandalism due to heavy recreational traffic.
Consequently, the highly portable, inconspicuous, inexpensive Gem was ideal for this campaign.
Preliminary results of this work are encouraging. Infrasound recordings are coherent across the network and a semblance grid search reveals a source region downstream of one of the hydraulic jumps ( Fig. 4.7). This project demonstrates the viability of dense networks of Gems for investigating and tracking fluvial hydrodynamics.

Do-It-Yourself Instrumentation
The rise of custom hardware has helped researchers in Earth science (Mallon, 2014) and other fields by reducing costs and increasing specialization of research equipment (Pearce, 2012). Custom instruments based on the popular Arduino platform are used for a range of tasks including detecting radiation, measuring pH, and analyzing DNA.
This trend has facilitated the sharing of open-source hardware, enabling researchers and labs to benefit from designs created by others while sharing their own designs in Figure 4.7: High-semblance source regions identified by grid search at Big Falls, South Fork of the Payette River, Idaho (Ronan et al., 2017). Infrasound data were band-pass filtered from 3 to 6 Hz before calculating semblance. The ability to image infrasound source regions in complex rapids using local instrument networks demonstrates this technique's usefulness for studying and tracking hydrodynamics in river features.
turn. New designs and ideas can spread more rapidly and cheaply than in a typical commercial marketplace.
Scientists who design their own instruments benefit from a thriving DIY community, whose members range from casual hobbyists to engineers. Several companies serve this community by providing easy-to-use components, detailed instructions, and online support. DIY-oriented suppliers used in this project include Adafruit, Evil- quickly. This free expert support is invaluable to beginning instrument designers.
Although DIY instrument design is more feasible now than in the past, scientists must also consider the disadvantages of designing their own instruments. The main disadvantage is the amount of time that must be invested in development and the engineering skills that must be mastered. Additionally, by building their own instruments, scientists take full responsibility for ensuring that instruments work once in the field. Rigorous instrument testing is a long process and detecting and correcting every possible malfunction is difficult, so scientists may prefer to leave such tasks to instrument engineers. Another potential difficulty faced by all non-standard or adjustable instrumentation (not just DIY instruments) is the need to conform to standard formats and completely account for instrument response when archiving data.
(Gem data files can be converted to common seismic formats easily, but researchers must keep track of any adjustments they make that affect the sensitivity or frequency response.) Finally, the manufacturing process is affected by economies of scale: components or manufacturing processes can be much more expensive or difficult when building only a few units (typical for an individual) than when building hundreds or thousands of units (more typical for an instrument company). All of these difficulties had to be overcome in the process of designing and building the Gem; however, in our case, we believe the outcomes of the project justify the challenges we faced.

Conclusions
We designed and built the Gem infrasound logger and collected useful data with it in the field. This instrument is considerably cheaper and more convenient for infrasound campaigns than more expensive commercial systems. Instrument design does not necessarily require the resources of an electronics company: the do-it-yourself electronics market offers easy-to-use components and free design software accessible to non-engineers. Designing and building their own equipment has helped scientists in a range of fields, and geophysicists can benefit from adopting this practice as well.
Future work will explore specialized applications of the Gem. We expect further innovations related to balloon-borne infrasound recording as this new data collection method matures. We also expect to see the development of novel means of infrasound recording. For example, we have begun testing the Gem as a payload for a multirotor UAV, and projects to record volcano infrasound by quadcopter are being planned. Additionally, we are exploring the possibility of recording infrasound underwater using the Gem. The ability to connect auxiliary sensors to the Gem creates further opportunities for specializations: for example, a Gem recording volcano infrasound could simultaneously record radiance of the vent, or a balloon-borne Gem could record absolute atmospheric pressure alongside infrasound. Finally, the Gem will continue to be used in traditional ground-based infrasound projects.

ERUPTIONS Summary
Volcanic ejecta contain a wealth of information on explosion and plume physics, but inferring source processes requires atmospheric effects to be understood and accounted for. This section introduces tephra transport modeling through the ambient atmosphere and discusses its application in two very different eruptions: the 14 July 2013 vulcanian eruption at Volcan Tungurahua (Ecuador) and the 3 March 2015 lava fountain at Volcan Villarrica (Chile). In the Tungurahua study, I show that near-vent acoustic sensors recorded infrasound from falling blocks (whose impacts were independently observed at least 2 km from the vent) and that the impact times and locations are consistent with reasonable eruptive parameters. In the Villarrica project, I show that a narrow (∼ 10 − 15 km wide, depending on downwind distance), elongate (at least 30 km) deposit of porous scoria lapilli can be explained with a tephra transport model that considers actual atmospheric properties from the time of the eruption, and that wind shear affects the deposit in subtle but noticeable ways. Although tephra properties vary widely and the physics of tephra transport is highly scale-dependent, a simple Lagrangian (particle-tracking) transport model is found to be an efficient and flexible way of modeling tephra trajectories in a variety of conditions. A cross-sectional area of particle (m 2 ) C D drag coefficient (unitless) g acceleration from effective gravity (true gravity plus centrifugal force) (m · s −2 ) G total gravitational, centribugal, and buoyant forcê k vertical unit vector, defined as parallel to effective gravity K diffusivity tensor (m 2 s −1 ) m mass of particle Q volumetric concentration of tephra (kg/m 3 ) r particle radius (m) r turbulent wind velocity vector (m · s −1 ) S source terms in conservation equations wind velocity vector (m · s −1 ) x cartesian position vector [x, y, z] (m) x, y, z cartesian coordinates (m) µ dynamic viscosity (P a · s) ρ density (kg · m −3 ) τ scale time of turbulent eddies (s) φ flux of a conserved quantity Q (units of Q times m −2 s −1 ) ∇ the 'del' operator ( ∂ ∂x , ∂ ∂y , ∂ ∂z )

Introduction
Pyroclastic fall deposits are tephra (volcanic ejecta) emplaced by falling to the ground through air after having been ejected explosively from a vent or carried upward in a plume (contrasting with volcanic deposits transported along the ground as lava, rockfalls, or pyroclastic density currents). They have been a major research topic in volcanology due to their ability to reveal properties of past eruptions (Houghton et al., 1999) and their role as a major hazard (Wilson et al., 2012).

Physics of Tephra Fall
Change in ejecta momentum (m ∂v ∂t ) of a particle depends on the combined gravitationalbuoyant forces acting on it (G, determined by particle size and mass), and atmospheric drag (D, determined by steady winds and random turbulent motion-discussed in section 1.3.1): Ejecta trajectories (with velocity v and position x) then depend on the initial position x 0 , initial velocity v 0 , and cumulative accelerations during flight:

Gravitational Force
Gravitational force (−mgk) on a particle is approximately constant in the atmosphere.
The buoyancy force (V tephra ρ atm gk) varies with atmospheric density depending on height, but is generally small due to the large density difference between air and rock.

Drag Force
Drag force depends on the particle and wind velocities as well as atmospheric properties. A particle is in equilibrium with zero acceleration when the drag force exactly balances the downward gravitational force; this occurs when the horizontal velocity of the particle is equal to the local horizontal wind speed and the vertical velocity of the particle is equal to its terminal velocity (or settling velocity) v s . Low-mass particles like ash typically reach equilibrium so quickly that their trajectories can be calculated accurately by following the horizontal wind velocity while falling at v s ; in this case, the initial velocity of the particle can be neglected (Folch, 2012). However, higher-mass particles like lapilli and blocks might take significant time to reach equilibrium or might never do so before impacting the surface, meaning that in such cases v 0 cannot be ignored.

Rayleigh Drag (Turbulent Flow)
The method used to calculate drag force depends on the dimensionless Reynolds number Re = ρatmvL µ , a ratio of inertial forces to viscous forces that determines whether flow is laminar (low Re) or turbulent (high Re) (Folch, 2012). In turbulent flow (high Re), viscosity is negligible. Rather, the drag force arises from Rayleigh drag (or form drag), the process of displacing air in the path of the projectile (Shuttleworth, 2012;Folch, 2012), defined as where C d is a dimensionless drag coefficient determined by the particle's geometry (0.6 for a sphere, 0.8 for a cube, and higher for rougher geometries) and A is the cross-sectional area of the particle. Rayleigh drag typically applies to coarse ash and larger particles (d > 100µm). A particle is said to be falling at its terminal velocity or settling speed when the upward drag force exactly balances the downward gravitational and buoyant forces mg(1−ρ atm /ρ tephra ); with Rayleigh drag, this balance happens at a settling speed of v s = 2mg C d ρatmA . Settling speed of blocks and bombs can be high enough that particles may never reach it before striking the ground; such particles typically follow ballistic-like trajectories (e.g., Mastin, 2001).

Stokes Drag (Creeping Flow)
On the other hand, particles moving through air with a low Reynolds number (Re < 0.1) are affected mainly by viscous forces, a process called Stokes drag (Slade, 1968): Where µ is the air's dynamic viscosity and r is particle radius. For very fine ash (d < 15µm), Stokes drag as defined in eq. (5.5) fails to account for effects of free molecular flow and the Cunningham slip correction must be applied (Slade, 1968).
Stokes drag mainly applies to fine particles (d < 100µm), which generally reach equilibrium rapidly due to their low inertia. Such particles fall at a settling speed of v s = 2gr 2 (ρ tephra −ρatm) 9µ .

Atmospheric Winds and Eddy Diffusivity
Atmospheric motions can generally be represented as large-scale ambient winds (vary- Instantaneous diffusive motions appear random due to their turbulent nature, and each particle's trajectory can be treated as Brownian motion following the Langevin equation (Boughton et al., 1987): where K and r i are the diffusivity and diffusive component of velocity (at time step i) in a given dimension, τ is the time scale of turbulent eddies (about 1 minute, as shown in figure 1.1), δt is the time interval, and N (0, 1) represents a normally distributed random number with mean zero and variance one.
For time scales longer than the characteristic turbulent time scale, turbulent motions can be represented well by deterministic diffusive motion following a Fickian flow law Qr = −K∇Q, where Q is the concentration of a tracer. The eddy diffusivity tensor K (discussed in section 1.3.1) is commonly written as a diagonal matrix in Cartesian coordinates with equal horizontal components and a different vertical diffusivity (often assumed to be zero) (e.g., Bonadonna et al., 2005b): Under certain assumptions (discussed in section 5.2.5), diffusive motion leads to an analytical solution where particle concentration follows a two-dimensional Gaussian distribution whose variance is σ 2 = 2K h t. (Folch, 2012)

Inertia-less Simplification
Due to the relatively low mass of ash and some lapilli, particles of those sizes are often modeled assuming the inertial term in eq. 5.1 to be negligible. In such a situation, the drag force must balance the gravitational-buoyancy force exactly. Because the gravitational-buoyancy force is purely vertical, this balance occurs when particles move at the same horizontal velocity as the surrounding air (so the horizontal drag is zero) and fall at a speed such that vertical drag force equals the gravitationalbuoyancy force (a speed called the settling speed or terminal fall speed):

Tephra Transport Modeling
Numerical models of tephra dispersal can follow two main approaches of solving (5.1).
"Lagrangian" or "particle-tracking" schemes track a set of particles from their origins to the ground by approximately solving (5.1) in discrete time steps. "Eulerian" models, on the other hand, discretize both time and space by dividing the atmospheric volume into a set of cells, and calculate the flow of a continuous tephra cloud from cell to cell. A third category of models, Gaussian models like Tephra2 (Courtland et al., 2012), reduce computational expense by using analytical solutions to the advectiondiffusion equation, but consequently are only valid under certain conditions.

Lagrangian Models
Lagrangian models predict particle trajectories by combining wind and settling velocities with turbulence as v = w − v sk + r (5.10) where v is particle velocity, w is wind velocity, v s is terminal fall speed, and r is random turbulent motion.
When a particle's inertia is non-negligible-such as bombs and blocks-a more complicated and general equation of motion must be used: where ρ atm and ρ tephra are the densities of the atmosphere and projectile, C d , A, and m are the drag coefficient, cross-sectional area, and mass of the projectile, and g is the acceleration from gravity.
Particle trajectories depend heavily on atmospheric properties including density, viscosity, and wind velocity; therefore, calculating tephra paths requires atmospheric structure to be measured, modeled, or assumed. Such atmospheric parameters are invalid very near the vent where hot, fast-moving gas erupts alongside the ballistics.
Failing to account for these effects has caused past studies to severely overestimate ballistic ejection speed and explosive pressure (Fagents & Wilson, 1993;Morrissey & Mastin, 1999), and it is necessary to account for the speed of eruptive gas if the purpose of the modeling requires use of the true initial velocity of particles at the vent itself. However, if the objective is simply to calculate realistic particle trajectories through ambient air and recovering the true initial velocity is not important, the role of erupting gases can be neglected because they affect only a small part of the particle's trajectory.
Lagrangian models track individual particles within a flow and therefore do not explicitly track the flow as a continuum. However, they have a distinct advantage in that they can be applied to particles with arbitrarily variable and complex physics with no simplifying assumptions.

Eulerian Models
Equation 5.1 can be reformulated into an Eulerian framework called the advectiondiffusion-sedimentation equation, enabling airborne tephra to be treated as a continuous fluid flow (e.g., the model FALL3D (Folch et al., 2009)). The advection-diffusionsedimentation equation can be derived starting with the fundamental conservation equation ∂Q ∂t + ∇ · φ = S, where Q is the conserved quantity (in this case, concentration of tephra particles), φ is the flux of the conserved quantity, and S represents sources or sinks of the quantity. Because flux of tephra occurs strictly due to its velocity (as opposed to terms like conduction or forces, which can appear in other conservation equations), φ = Qv. Eulerian models are mainly applied to particles with negligible inertia, meaning that particle velocity is determined by gravitational settling, wind, and eddy diffusion: φ = Qw − Qv sk − K∇Q. This leads to the advection-diffusion-sedimentation equation where Q is the tephra concentration, w is the wind velocity, v s is the terminal fall speed, K is the diffusivity tensor, and S is the source term (representing injection of tephra into the atmosphere by the volcano). For a heterogeneous set of tephra, this equation must be solved several times for different particle types because terminal settling velocity is highly variable among particles of different size, shape, or density.
Eulerian models are most useful when the flow being modeled consists of a large number of particles with similar behavior, and has the advantage of explicitly representing a flow as a continuum. Although it would be convenient to represent the flow of tephra as a continuum in the 2015 eruption of Villarrica, I use a Lagrangian model instead because of its flexibility in modeling particles with different properties.

Gaussian Models
An analytical solution to the advection-diffusion-sedimentation equation can be applied under certain conditions. This method makes the following assumptions: zero vertical wind or diffusion, horizontal winds varying with elevation only, constant diffusivity tensor, negligible divergence of settling velocity, and an instantaneous point-source release of particles. Under these conditions, the advection-diffusionsedimentation equation simplifies to where x 0 , t 0 , and S are the location, time, and magnitude of the particle source (Folch, 2012). This system has the analytical solution where the center-of-mass position varies as z c (t) = t t 0 v z (z c (τ ))dτ , x c (t) = t t 0 v x (z c (τ ))dτ , and y c (t) = t t 0 v y (z c (τ ))dτ . This analytical solution has the form of a two-dimensional Gaussian distribution whose mean is the center of mass and whose variance increases with diffusivity and time.

Methods
I used a custom tephra transport model (Anderson, 2018) written in R for ease of integration into data analysis scripts (for example, the code provided in Appendix 6).
The transport model is a Lagrangian (particle-tracking) scheme that releases particles at an initial velocity and position at the vent or in the plume, then tracks them as they fall through the atmosphere until they cross the topographic surface. The model makes no assumptions about particle inertia or diffusivity, and is therefore equally valid for tracking blocks following ballistic trajectories or windblown, diffusing coarse ash. However, it does assume that the Reynolds number is high enough for Rayleigh drag (eq. 5.4) to apply; this assumption is valid for the coarse tephra studied in this dissertation (Folch, 2012). I use atmospheric models obtained from the Global Forecast System Analysis dataset (GFS-ANL) (Kalnay et al., 1990) using the rNOMADS package to download data (Bowman & Lees, 2015a). The GFS data have a temporal resolution of six hours and spatial resolution of 0.5 degrees of latitude and longitude and include vertical profiles of pressure, temperature, elevation, and wind (expressed as a two-dimensional vector in the horizontal plane). Because the latitude and longitude of Volcán Villarrica do not fall exactly on an atmospheric model location, the atmospheric model used for tephra dispersal modeling is interpolated between the four nearest GFS-ANL atmospheric models. The vertical resolution of the model is 5000 Pa in pressure coordinates, corresponding to spacing of about 500 m and 1000 m at elevations of 1500 and 8000 m above sea level. I calculate atmospheric density using the provided pressure and temperature and the ideal gas law.
GFS atmospheric models are representative of the atmosphere as a whole but are not accurate immediately around the vent during an eruption, and may fail to capture topographic effects near the ground. I do not attempt to model tephra transport in the plume rising above the vent (instead assuming that tephra leaves the plume at elevations ranging from the vent to the observed top of the plume) and assume that topographic effects on near-surface winds have negligible effects on tephra.
Topography around these volcanoes was taken from digital elevation models (DEMs) in the ASTER dataset (NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team, 2009).

Tungurahua, Ecuador
Infrasound recordings from the vulcanian eruption at Tungurahua discussed in chapter 3 include an unusual phase following the main arrival, characterized by an excursion in source direction away from the vent and subsequent return to the vent.
Such a rapid source migration occurring so early after the eruption is difficult to explain with already-known acoustic source types at volcanoes. Instead, I attribute this phase to impacts of falling ballistics.

Explaining a rapidly moving volcanic acoustic source
Infrasound recordings from spatially separated microphones can be compared to determine time lags among the sites. The microphones at the Tungurahua stations were installed as linear arrays oriented toward the vent, meaning that analysis of time lags can determine the angle between the wave incidence direction and the vent.
Therefore, an angle of zero degrees (as recorded during the main blast wave arrival) corresponds to waves arriving directly from the vent, and an angle of 90 degrees corresponds to waves arriving from either side of the array or above it.
Changing incidence angles of infrasound recordings reveal a new phase following the main blast. In this phase, the angle at station 9025 deviates from 0 • (in-line with the vent) around 10-15 s, reaches 90 • at 19 s and peaks slightly later at 106 • , and again approaches 0 • around 24 s ( fig. 5.2a). The same signal also appears on the lower station A3H with a 7.6-s delay.
I consider four types of potential acoustic sources with a high incidence angle (relative to the vent) to explain this unusual signal: a plume, pyroclastic density currents (PDCs), echoes, and falling ballistics. Although volcanic plumes have been identified as acoustic sources at other volcanoes , I consider a plume to be an unlikely explanation for this particular because it would need to expand at an unreasonably fast speed (averaging at least 97m · s −1 ) in order to travel the 1850-m horizontal distance to the station in 19 s. Similarly, moving acoustic sources (e.g., rock falls and PDCs) are common at volcanoes and can be identified by their changing propagation directions (Ripepe et al., 2010;Johnson & Ronan, 2015).
However, an unreasonably fast speed of at least 120m · s −1 (about three times as fast as the fastest PDCs recently recorded at Tungurahua (Hall et al., 2013(Hall et al., , 2015) would be needed to reach the array in 19 s. Finally, I consider topographic echoes to be an unlikely source because no major reflectors have an appropriate two-way distance (6.5 km to explain an arrival at 19 s at a sound speed of 340m · s −1 ).
Bombs and impact craters (figure 5.1) were observed in the vicinity of the upper station and between the two stations after the eruption (Ortiz, 2013), motivating the investigation of falling ballistics as a possible acoustics source. Having discarded other explanations for this sound, I therefore test whether ballistic impacts (a source not previously documented in volcano infrasound studies) can explain this signal, and find good agreement between ballistic models and infrasound observations. I ran several explosion simulations, each testing a different initial speed and particle size combination. Each simulation tested a range of ejection angles, but initial speed and other ballistic properties were identical for all particles within a particular simulation. Particle radius r and initial speed v 0 varied between simulations, but for simplicity, I fixed particle density ρ tephra at 2500kg/m 3 and drag coefficient C d at 0.6 (greater than a sphere, but less than a cube), and assumed particles to be approximately spherical so that cross-sectional area A ≈ πr 2 . I performed several such model runs with different radii and initial speeds in search of a good fit to infrasound propagation data.

Ballistic Modeling at Tungurahua
Among the scenarios tested, several resulted in bombs falling first uphill of the array, then passing downhill of it, then moving uphill of it again ( fig. 5.2b). For Ballistic fall is the only candidate source mechanism that can explain this infrasound, and it requires no unlikely assumptions to do so. Although ballistic-derived infrasound had not been identified before this study, it is probably common in vulcanian eruptions but difficult to resolve with distant instrumentation; and future near-vent array studies may detect it as well.

Villarrica, Chile
After weeks of strombolian activity, Villarrica erupted in a lava fountain in the early morning of 3 March 2015. The lava fountain reached 1500 m above the summit and produced a plume rising to 6-8 km above the summit Romero et al., 2018;Sennert, 2015). Porous scoria lapilli coated the summit and eastern slope of the volcano ( fig. 5.3) and fell in a narrow band extending at least 20 km toward the east and southeast.

Field and Satellite Observations of Tephra
The dimensions of the tephra deposit were constrained through a combination of field mapping and remote sensing. Tephra from the March 2015 eruption could be distinguished from older deposits in the field by the glossy black appearance of fresh   scoria, and in many sites the presence of intercepted fresh lapilli on tree branches confirmed the difference ( figure 5.4). I mapped the presence or absence of tephra at several field sites in order to constrain the width of the deposit. Additionally, I recorded the dimensions of a sample of at least ten largest pieces found at certain sites; these sample dimensions were recorded either as scaled photos (figure 5.5a) or as three-axis dimensions recorded in the field.
Unfortunately, the ability to completely map the tephra deposit was limited by access to field sites, so I complemented my field observations with satellite images before and after the eruption (figure 5.6). In many locations covered with snow or ice, presence or absence of tephra could be determined easily in satellite imagery. My resulting tephra fall map shows a generally narrow deposit extending at least 30 km east to southeast from the vent (figure 5.5b).
As expected, particle size grades in the downwind direction: sites far from the vent have coarser tephra than sites near the vent. Additionally, grading is seen in the cross-wind direction: at two sites at approximately equal distances from the vent (12-13 km), tephra was significantly coarser along the deposit's southern boundary than along the northern boundary (figure 5.7). I did not observe site-dependent variation in tephra characteristics other than size among the sites visited.

Modeling Tephra Dispersal
I modeled the Villarrica plume and deposit using the custom transport model de- Two parameters that vary among tephra particles should determine their expected    Wind speed and direction varies with elevation; most importantly, winds below 6000 m are slower and around to 90 − 110 • , while winds higher than 6000 m blow faster and around 120−130 • . Consequently, tephra that falls from high in the plume will be advected more to the southeast, while tephra that falls low in the plume will be advected more to the east.
trajectory: the height at which they leave the eruptive column (and begin drifting with the wind and falling at the settling speed), and tephra aerodynamic parameters like cross-sectional area, drag coefficient, and mass. I test a range of plausible values for both of these two properties, without attempting to estimate their distribution.
The 3 March 2015 plume at Villarrica was estimated to reach an elevation of 6-8 km above the crater (itself at 2850 m above sea level) (Sennert, 2015). Consequently, I test particle initial elevations at 500-m increments from 3000 m (approximately the elevation of Villarrica's vent) to 11500 m above sea level. All particle trajectories begin above the vent at their initial elevation, moving at the same velocity as the wind.
The aerodynamic properties of the tephra are generally unconstrained, mainly because the drag coefficient and frontal cross-sectional area are difficult to measure in deposit due to their highly irregular shapes (e.g., figure 5.4). Rather than attempt to estimate these values, I instead consolidate them into a single "ballistic coefficient" the ballistic coefficient is a measure of a particle's resistance to Rayleigh drag. Due to lack of rigorous constraints, I consider ballistic coefficient to be a somewhat free parameter, though it must vary substantially due to the range of tephra sizes found in the field. I found that a ballistic coefficient range between 0.42 and 4.2 kg/m 2 resulted in particles falling in distances relevant to the tephra deposit observations. This range is equivalent to smooth spheres (C d = 0.47) of dense rock (2700kg/m 3 ) with particle diameter ranging from 0.11 to 11 mm. Considering that the Villarrica tephra had considerably lower density (due to its high porosity), higher drag coefficient (due to its irregular shape), and probably higher ratio of frontal area to volume, I expect this range of ballistic coefficients to correspond to much coarser particles, consistent with my observations in the field. I first calculate particle trajectories in the absence of turbulence (figure 5.9). Without turbulent diffusion, these trajectories are completely deterministic: all identical particles released at equal elevation will follow the same trajectories. Although the real atmosphere does have turbulence, which significantly affects the width of the deposit, it is illustrative to test deterministic trajectories first.
In a homogeneous wind field without turbulence, tephra impact distance would be equal to the wind velocity times the fall time of the particle (itself dependent on the initial elevation of the particle and the ballistic coefficient). The resulting deposit would be a thin straight band with particle size decreasing with increasing distance. However, the deposit in figure 5.9 varies in both orientation and width. This can be explained by the variability of the wind direction with elevation: at close distances to the vent, the deposit aligns with the wind direction at low elevations (about 100 • ), while at greater distances the deposit aligns with the high-level winds (about 120 • ). Intermediate distances (around 5-15 km) include fine particles that fell from low elevations and coarse particles that fell from high elevations; consequently, both trends are present here, making the deposit wider and graded south-to-north. Next, I model the tephra fall including the diffusive effects of atmospheric turbulence. Because the tephra remained in the atmosphere for a long time, atmospheric turbulence should also cause variation in trajectory for identical particles that leave the plume at the same height. I use the turbulent diffusivity tensor K to represent turbulence in both the plume and the ambient atmosphere (eq. 5.12). Because K is poorly constrained a priori, I consider it a free variable and test different values in order to match the observed deposit width. I find that a horizontal diffusivity of 750 m 2 /s can reproduce the shape and width of the deposit well (figures 5.10, 5.11). This diffusivity value is within the normal range used for atmospheric diffusion modeling (∼ 10 1 −10 4 m 2 /s), including modeling volcanic plumes (Bonadonna et al., 2005a). Effects of wind direction shear are again evident in this model.

Discussion of Villarrica Tephra Modeling
Early observations of the 3 March 2015 tephra deposit at Volcan Villarrica naturally prompted questions on why it took a long, narrow form. In this section, I have demonstrated that the deposit can be reproduced by simple atmospheric transport modeling that includes the atmospheric conditions at the time of the eruption, the local topography, the observed plume height, and reasonable values of ballistic coef- Figure 5.9: Modeled tephra deposit assuming no turbulence. Point size and color correspond to ballistic coefficient: small, blue points represent particles with ballistic coefficients (generally finer tephra) and large, red points represent large ballistic coefficients (generally coarser tephra). Particles grade downwind as expected. Additionally, multiple effects of wind shear become evident 5-15 km east of the vent: around that distance, the axis of the deposit changes direction, the width increases, and cross-wind grading is evident. Figure 5.10: Modeled tephra deposit assuming that atmospheric turbulence results in an eddy diffusivity of 750 m 2 /s). Point size and color again correspond to ballistic coefficient: small, blue points represent particles with small ballistic coefficients (generally finer tephra) and large, red points represent large ballistic coefficients (generally coarser tephra). Due to diffusive effects of atmospheric turbulence, the deposit is noticeably wider than in the turbulence free model shown in figure 5.9. However, similar effects of wind shear remain visible here: the axis of the deposit varies and ballistic coefficient grades from south to north. Figure 5.11: Snapshots of tephra transport models at Villarrica. Left column assumes no diffusion (particles advected by steady winds only); right column assumes horizontal eddy diffusivity of 750 m 2 s −1 . Particles disperse more widely in the models including diffusion by turbulent eddies. However, tephra grading downwind to the southeast (due to large particles settling faster) and cross-wind to the northeast (due to wind shear) is evident in both models. ficient and eddy diffusivity. Additionally, I have shown that wind shear played an important role in the transport process: a small but noticeable shear in wind direction (20 − 30 • ) between about 5000 and 8000 m above sea level can explain particle size grading in the cross-wind direction and may have caused variability in deposit width and azimuth. It was not necessary to consider complex and poorly constrained (but potentially important) processes like plume convection, particle-particle interactions, particles breaking or aggregating within flight, or topographic effects on boundarylayer winds.
Limitations of this work include uncertainty related to limited time and access for field mapping, the unclear relationship between particle size and ballistic coefficient, and the lack of constraints on eddy diffusivity along particle trajectories. Although the field map of tephra presence and absence constrains deposit width and orientation, significant parts of the deposit were simply inaccessible, and it was not possible to systematically record tephra mass load or particle sizes at all sites. Further, the difficulty in determining aerodynamic properties of tephra (due to their irregular shape, uncertain orientation, and possibility of having broken during flight) prevents us from inferring particle ballistic coefficients at sites based on tephra morphology.
Finally, eddy diffusivity remains a free parameter within a wide range of reasonable atmospheric values (10 1 − 10 4 m 2 /s). If more constraints were available on eddy diffusivity and ballistic coefficient, particle trajectories could be modeled with greater certainty, enabling inferences to be made about plume processes and the distributions of tephra sizes and initial elevations within the plume.

Conclusions
Pyroclastic fall and ballistics matter to the volcanology community due to the information they carry about eruptive processes and their roles as hazards. Infrasound data recorded in the 14 July 2013 eruption of Volcan Tungurahua (Ecuador) show that near-vent infrasound is a viable means of detecting falling ballistics without depending on visibility or placing human observers in harm's way. Ballistic trajectories are useful because they can reveal explosive processes like initial gas pressure. Future work could use comprehensive telemetered near-vent infrasound installations to identify impact times and locations in order to recover explosive parameters in near real time (tens of seconds).
Pyroclastic fall in the 3 March 2015 eruption of Volcan Villarrica (Chile) fell in a long, narrow band extending from the vent to the east and southeast. Transport modeling based on the actual atmospheric properties at the time of the eruption succeeded in reproducing the general shape of the deposit. Additionally, wind direction shear was found to cause particle size grading in the cross-wind direction.
Although many valid tephra transport models already exist, I found it advantageous to use a custom code (Anderson, 2018) to model tephra trajectories. Most importantly, using software written in a common scientific programming language (R) enabled me to easily and flexibly integrate it into scripts that facilitated visualizing results and testing a range of eruptive and atmospheric conditions. Additionally, using an easily modifiable code enabled me to extend the model to consider additional features beyond the original purpose; consequently, the model is equally valid for calculating trajectories and flight times for ballistic blocks as for fine lapilli transported by wind, settling velocity, and eddy diffusion. The resulting code reproduced the ob-served tephra deposit at Villarrica and would be useful for modeling future eruptions, with the potential for real-time hazard forecasting applications. Shock waves occur when the initial condition contains strong pressure contrasts or supersonic flow (chapter 2). Shock wave behavior is complicated: celerity varies along the wave causing it to change shape while propagating, energy is lost from the wave and deposited in the air behind it, and the shock speed decreases as the shock decays. No simplifications of the Euler equations can be made when modeling shocks, and no analytical solutions are available. Numerical modeling of the full Euler equations requires complex and computationally expensive finite-volume methods.
Powerful volcanic eruptions produce shock waves, and interpreting them is difficult due to their nonlinearity, lack of analytical methods, and expense of modeling; in the absence of linear inversion methods, source inferences would require an expensive iterative inversion with a new simulation in each step. However, due to the scaling properties of the Euler equations, a single solution can be adapted to several different eruptive scenarios by application of scaling laws (eq. 2.8-2.9). This combination of numerical modeling and scaling laws provides a framework by which source parameters of explosive eruptions (in particular, the energy of the explosion) can be estimated at reasonable computational expense.
In problems involving infinitesimal pressure and velocity disturbances over length scales up to several km, the Euler equations are dominated by the pressure gradient and particle acceleration, resulting in acoustic waves (called infrasound at frequencies below 20 Hz) following the linear acoustic wave equation (eq. 1.32). Acoustic waves have simple physics: their speed and attenuation are independent of amplitude, their propagation is linear can be represented by convolutions, and they can be implemented in simple finite-difference models. Consequently, infrasound has a variety of powerful analytical and numerical methods that are useful for studying volcanoes. In chapter 3, I use a numerical model of the acoustic equations to analyze recordings of a powerful vulcanian explosion at Volcan Tungurahua, Ecuador. Results show that total erupted volume (gas and tephra) was around 0.5 km 3 , with excellent agreement between microphones located over 1 km apart. Additionally, infrasound signals reveal a short (0.7-s) displacement of air immediately preceding the explosion (interpreted as an uplift of the vent surface), as well as a complex period of tremor following the explosion including three types of vent activity. Placing sensors near enough the vent enabled infrasound to be recorded with relatively little topographic scattering, enabling detailed inferences of vent activity to be made.
Field data collection is an important part of volcano infrasound study, and requires instrumentation appropriate for the task. Chapter 4 describes the development, char-acteristics, and use of a new infrasound instrument designed for low cost, low power consumption, and portability. Instrumentation is an essential component of fieldwork and partly determines the kind of work that can be done and the kind of data that can be collected; instrumentation whose characteristics are not well aligned with the purpose of fieldwork can limit the success of field campaigns. I was motivated to develop a new infrasound logger because the price, weight, and power needs of existing data loggers limited the number of infrasound sensors I could install in volcano fieldwork, and restricted the installation sites in which I could place them. The resulting Gem infrasound logger can be purchased and transported in large numbers due to its low cost, weight, and power consumption, and can be installed in arbitrary geometries.
It has proven useful in volcano fieldwork as well as other applications; in particular, its light weight makes it an ideal payload for low-cost stratospheric balloons.
Finally, atmospheric transport of tephra is also governed by the Euler equations, but the result has very different behavior from pressure waves. Chapter 5 describes Lagrangian modeling of tephra from a lava fountain eruption at Volcan Villarrica, Chile. In this problem, the time scales are much greater than time scales relevant to atmospheric turbulent eddies (around 1 minute), meaning that turbulence can be treated as diffusive random motion. Additionally, the duration of plume transport in the area of interest is much less than synoptic time scales (hours to days), meaning that synoptic scale winds controlled by the pressure gradient and Coriolis force dominate and are approximately steady. Consequently, tephra transport consists of two simultaneous behaviors: advection by ambient winds and settling speed, and diffusion due to random motion from turbulent eddies. My Lagrangian model of tephra transport accurately reproduces the observed tephra deposit from the eruption, which consists of a long, narrow band with coarser tephra on the south side of the deposit.
I applied the same model to predict ballistic trajectories from the vulcanian eruption at volcan Tungurahua, Ecuador, finding that ballistic impact times and locations were consistent with unusual infrasound recorded at that time (which could not be explained using typical volcano infrasound sources).