Modeling source parameters of quasi-periodic tremor

Modeling harmonic tremor recorded at volcanoes is an essential practice in improving eruption forecasting methods and warning systems. We model the conduit dynamics of quasi-periodic tremor (chugging) recorded at Karymsky, Tungurahua, and Fuego volcanoes to estimate its source characteristics. Chugging mechanisms are estimated using two theoretical models originally derived in Garcés [1997, doi: https://doi.org/10.1029/1999JB900096] and Girona et al. [2019, doi: https://doi.org/10.1029/2019JB017482], respectively. Comparisons of the conduit and fluid output parameters suggests that chugging is primarily limited to near-surface oscillations and outgassing due to an accumulation of gas between eruptive episodes. The modeled results indicate clustered release of volatiles exsolved from a deeper magma conduit region, triggered by an initial explosion. This interpretation is consistent with both infrasonic and seismic observations at each volcano.


Introduction
The detection and modeling of harmonic tremor is an important part of understanding volcano fluid dynamics and improving source process prediction methods [Roman 2017]. Quasi-periodic tremor, otherwise referred to as chugging, is a type of harmonic tremor generally characterized by a discrete series of emergent acoustic pulses following in the coda of a highfrequency initial explosion (approximately 10-50 s time lag). Chugging signals have a common dominant frequency ranging 0.6-1.4 Hz accompanied by harmonic frequencies not typically observed in similar jetting events or long-duration explosions. A single chugging sequence can last up to 90 s and in some cases behaves intermittently with several (up to 4) discrete emergent sequences occurring in the coda of a single explosion . Figure 1 shows a characteristic chugging sequence, as well as the accompanying initial explosion, for an episode recorded at Karymsky volcano, Russia, in August, 1997. Both the acoustic (top) and seismic (bottom) data are shown to demonstrate that some characteristics of chugging, such as double-periodicity and asymmetry, are better resolved in the infrasound signal. Thus, only infrasound data are modeled and discussed in this paper.
Acoustic chugging has been observed and recorded at several volcanoes around the world including Karymsky, Russia , Tugurahua, Sangay, and Reventador, Ecuador [Johnson and Lees 2000;Lees and Ruiz 2008;Ruiz et al. 2006], Arenal, Costa Rica [Benoit and McNutt 1997] and Fuego, Guatemala [Lyons et al. 2009]. Despite many geological and compositional Corresponding author: juliant@live.unc.edu differences among these volcanoes, chugging events recorded at each site share nearly all of the frequency and physical properties unique to this type of tremor.
We model chugging data from three volcanoes-Karymsky, Tungurahua, and Fuego-using methods developed by Girona et al. [2019] and Garcés [1997] to estimate the appropriate source parameters and mechanisms able to resolve the features of chugging described in this section. Accurate modeling of harmonic tremor is crucial to the development of new eruption forecasting methods and creating better warning systems. Chugging at Fuego volcano, for example, has been linked to precursory activity indicative of paroxysmal, long-duration eruptions lasting 24-48 hr Lyons et al. [2009]. Volcanic tremor at all three volcanoes has been reported to precede nearly all recent eruptive periods during which warning systems and evacuation plans were executed [Johnson and Lees 2000;Lyons et al. 2009;Ruiz et al. 2006]. The models used in this research are highly variable in geometry and source process, thus providing a basis for comparison of the physical attributes of simulated results and an effective method for discerning an interpretation of chugging source dynamics. There are currently several models that have been developed and used to estimate and predict the source mechanisms driving the spectral and temporal signatures of harmonic tremor recorded at many active volcanoes around the world. Crack extensions and deep conduit hydraulic pressure due to magma flow have been proposed to describe the source of harmonics observed in acoustic or seismic waves [Aki et al. 1977;Julian 1994]. These models suggest that resonance within a deep elastic channel occurs due to the fluid-elastic interactions of the magma with the Modeling source parameters of quasi-periodic tremor Traphagan & Lees, 2020 fluid boundaries; thus, nonlinear particle oscillations are suggested to be excited into motion by the elasticity of the channel walls. More recent studies [Johnson and Lees 2000;Lees and Bolton 1998;Lees et al. 2004] have suggested that harmonic tremor, especially chugging episodes, may instead be driven by surface degassing or shallow conduit processes, instead of deep fluid flow dynamics as predicted by some models. This research focuses on the simulated model output of the theoretical models proposed in Girona et al. [2019] and Garcés [1997] in order to consider the differences between a conical conduit body and near-surface gas oscillation model and to comparatively interpret the model results for an appropriate chugging source model. A comprehensive model would be able to successfully resolve the spectral characteristics of chugging while remaining consistent with measurements and direct observations of the physical properties of the volcano.

Data and Study Regions
Acoustic quasi-periodic tremor signals recorded at Karymsky volcano, Russia, Tungurahua volcano, Ecuador, and Fuego volcano, Guatemala are used to model chugging signals as a method of constraining their source parameters. These three stratovolcanoes are similar in composition-viscous andesitic to basaltic magma with high gas content-and exhibit similar eruptive patterns, which are typically Strombolian in type with frequent continuous degassing [Braitseva 1990;Hall et al. 1999;Johnson et al. 2003;Martin and Rose 1981;Ozerov et al. 2003;Roggensack 2001]. These sites were chosen for this study due to frequent visual, audio, and recorded observations of chugging episodes at each volcano [Johnson et al. 1998;Johnson and Lees 2000;Lees et al. 2004;Lyons et al. 2009;Ruiz et al. 2004;Steele et al. 2014]. Maps of the volcanic regions are provided in Figure 2.

Karymsky
Recent eruptive activity at Karymsky is characterized mostly by small explosions and degassing. Karymksy is a stratovolcano surrounded by lava flows typically less than 200 years old, as well as a deep hydrothermal fluid system positioned along the fault patterns and intersections of the basement rock [Leonov and Grib 2005].
Eruptions of Karymsky are Vulcanian or Strombolian in type with plumes rising up to 17 km and typically consist of volcanic bombs, ash, and viscous andesitic lava flows [Braitseva 1990;Ozerov et al. 2003]. The cone of the volcano is located in a system of several overlapping calderas ranging to mid-Pleistocene in age [Erlich 1986]. Tephra from the youngest cone has been carbon dated at about 5.3 ka [Braitseva 1990], which is situated within a caldera about 7.7 ka in age and 2.5 km in radius [Hrenov et al. 1982]. Karymsky is the most active volcano in eastern Kamchatka Peninsula along the region's coastal volcanic belt. The volcano has undergone two major eruptive cycles separated by about 2300 years. The most recent volcanic activity began in 2001 and continued steadily until 2016 with sporadic episodes until September 2016, whereas the prior eruptive episode spanned nearly 12 years between 1970 and 1982 followed by a brief 4-year active period starting in January 1996 [Ozerov et al. 2003].

Tungurahua
Tungurahua is one of Ecuador's most active stratovolcanoes and is located in the Ecuadorean Andes about 140 km south of the capital city of Quito. Tungurahua's cone has steep flanks of about 30°and an open crater at the summit [Ruiz et al. 2006]. Eruptions are Strombolian in type and produce andesitic and dacitic magmas that form into viscous lava flows and pyroclastic flows. Recent eruptive activity at Tungurahua has produced plumes over 7 km [Ruiz et al. 2006] and pyroclastic flows killing at least five people in neighboring towns. Continuous degassing and jetting signals have also been recorded during seismoacoustic deployments since 1999 [Johnson et al. 2003]. Historical eruptive activity at Tungurahua includes five major eruptive periods since the 16th century, the most recent of which began in 1999 and continued through September 2016 [Hall et al. 1999]. The structure of the volcano is comprised of three edifices dating 14 ka for the oldest and nearly 3 ka for the most recent caldera collapse. The present cone has built to about 50 % of the previous cone size before its collapse and has mostly generated pyroclastic eruptive products in the last thousand years.

Fuego
Fuego volcano is a stratovolcano located in central Guatemala with frequent and historically sub-Plinian and Strombolian eruptions. The Strombolian acitivity at Fuego typically occurs over long durations of months or, in some cases, years. Most recently, eruptive activity began in 1974 as a series of four intense (VEI 4) sub-Plinian explosions over a period of a month followed by prolonged low intensity Strombolian activity that lasted into 1979 [Martin and Rose 1981]. Since Figure 2: Topographic maps of Karymsky (top), Tungurahua (middle), and Fuego (bottom) volcanoes with station locations. Contour intervals are 10 m for Karymsky and 50 m for Tungurahua and Fuego.

Presses universitaires de �rasbourg
Page 253 Modeling source parameters of quasi-periodic tremor Traphagan & Lees, 2020 then, Fuego's eruptive activity has been characterized by smaller (VEI 1-2) sub-Plinian eruptions occurring between 1980 and 2000, and continued activity from 2002 to the present. Currently, the volcano has produced frequent lava and pyroclastic flows as well as lahars and near-daily small scale eruptions [Lyons et al. 2009]. Lava flows from Fuego are typically basaltic or basaltic-andesitic in composition and contain high volatile contents of mostly H 2 O [Roggensack 2001].

Seismoacoustic Deployments
Acoustic signals were recorded at Tungurahua volcano by a network of 5 co-located broadband seismometers and infrasound microphones for 8 days in May of 2010. Seven co-located seismic and acoustic sensors were deployed around the vent of Fuego volcano, as well as an acoustic array of 4 microphones grouped approximately 23 km northeast from station VF03 (see Fig

Modeling Methods
The methods derived in Girona et al. [2019] for an oscillating pocket of trapped gas beneath a permeable plug are used to model the parameters of chugging signals recorded at Karymsky, Tungurahua, and Fuego with shallow source processes. We also model these chugging data with equations derived in Garcés [1997] and Garcés et al. [2000] as a comparative method for estimating harmonic properties of chugging signals represented as longitudinal modes through a conical conduit body.

Layered Conduit
Resonance in the magma conduit is described in Garcés [1997] as longitudinal wave propagation through a three-layered medium, each with a unique set of physical properties (i.e. viscosity, magma density, sound speed) characteristic of the flow regime of the conduit section. The deepest conduit section (S 3 ) is described as a region of homogeneous, single-phase flow with constant fluid injection due to decompression of a deeper fluid source [Garcés 1997;Garcés and Hansen 1998]. The flow regime in the intermediate layer (S 2 ) is characterized by a two-phase fluid composed of rising magma from the basal layer and volatiles exsolved from at the S 3 -S 2 boundary. The geometry of the conduit is ostensibly conical, an approximation that allows for the wave propagation to be described using Helmholtz resonance theory; the S 3 Figure 3: Simplified diagram demonstrating the geometry and mechanism for the three-sectioned conduit resonance model Garcés [1997]. Geometric variables include the conduit section length L i , cross-sectional area S i , and distance from the vent to the recording microphone r. Fluid properties such as the magma density ρ i , viscosity µ i , and sound speed c i are unique to the flow regimes in each section of the conduit. Figure adapted from Figure 1 of Garcés [1997].
layer is the widest, decreasing in width approaching the uppermost S 1 section. At depth z " 0, the resonating waves in the S 1 section are released into the atmosphere through a rectangular vent with dimensions aˆb and propagated to a receiver [Papale and Dobran 1994]. The equation for pressure recorded at microphone r meters from the vent is given by: where j " ?´1 , k 0 " ω{c 0 , ρ 0 " 1 kg m´3, and c 0 " 320 m s´1 are the wavenumber, density and sound speed of the atmosphere at the vent respectively. f ω is an angle-distribution factor estimating the amplitude of the acoustic signal propagating away from the center of a rectangular vent with dimensions a and b in a direction defined by θ "´z{r for all z [Garcés 1997;Morse and Ingard 1968]. Resonance in the conduit is driven by transient fluid injection from below the S 3 layer with velocity: where ν is the maximum source particle velocity, is Presses universitaires de �rasbourg

Page 254
Volcanica 3(2): 251 -262. doi: 1 .3 9 9/vol. 3. 2.251262 the exponential decay constant, and n is an integer. Fluid velocity throughout the conduit is determined by: where Z ij is the acoustic impedance contrast at the boundary of contiguous layers i and j, ρ i is the section magma density, c i is the sound speed, and the modal structure (M) of the magma conduit is a function of the the impedance conditions and attenuation coefficients for each conduit layer [Garcés 1997]. Thus, M defines the resonance structure of each section of the conduit, incorporating the effects of viscous attenuation via a complex wavenumber: where the absorption coefficient: is a function of the molecular relaxation time, τ i [Garcés 1997;Medwin and Clay 1998].
The attenuation and wavenumber are calculated over all three sections of the conduit across all frequencies to determine frequency and intensity of the resonant modes. The source particle velocity is used to establish the amplitude envelope of the synthetic signal, which is proportional to the transient mass flux input in the magma reservoir [Garcés et al. 2000].

Plugged Gas Pocket
The shallow source model derived in Girona et al. [2019] is used for comparison to the three-layered conduit model from Garcés [1997]. Oscillations are generated within a trapped pocket of accumulated gas beneath a permeable rock layer and exsolved from a deeper multiphase fluid source Girona et al. [2019]. Rising packages of bubbles burst upon impact on the gas pocket, driving the pocket into resonance. The model simulates each burst as a wavelet imposed on the synthetic time series output corresponding in time to the instant of impact and in amplitude to the bubble mass. The dominant oscillation frequency of the gas pocket is represented by: where the damped harmonic oscillator coefficients and constants (Γ 1,2 and γ 0,1 ) are auxiliary parameters related to the dimensions and pressure evolution of the gas pocket [Girona et al. 2019]. Pressure oscillations are calculated non-dimensionally using: where P ex is the external pressure, P 0 is the steady-state pressure, q m is the bubble mass, t m is the instant of bubble impact, N is the number of gas impulses, Γå nd ω˚are necessary auxiliary parameters for Taylor Series expansion, and ∆t m " t´t m and iterated using the Heaviside function Hp∆t m q [Girona et al. 2019]. The functions r c and r s represent hyperbolic cosine and sine functions when 4Γ 0 Γ 2 ă Γ 2 1 . The functions instead become cosine and sine functions when 4Γ 0 Γ 2 ě Γ 2 1 .
The resonance wavelet for the individual mass impulses is a function of the angular frequency (ω) and is represented by: which describes the frequency, amplitude, and envelope of the oscillations observed in the model waveform. To develop an amplitude envelope that resembles chugging, bubbles are simulated to impact the pocket at the instant of each chugging pulse (peak),

Presses universitaires de �rasbourg
Page 255 Modeling source parameters of quasi-periodic tremor Traphagan & Lees, 2020  [Girona et al. 2019] and the bottom row shows the spectra calculated by the resonating layered conduit model [Garcés 1997]. Columns are organized and labeled by volcano and red dashed lines indicate the dominant frequency of the recorded signal.
with mass proportional to the amplitude of the pulse. Impact times are determined based on an amplitude threshold (0.001 Pa) below which fluctuations in the data can be attributed to noise.

Numerical Parameter Optimization
We apply a root-mean-squared residual (RMS) condition to a constrained Nelder-Mead (NM) optimization algorithm to determine numerical estimates of tremor source parameters. This optimization method was chosen to numerically minimize these volcanic tremor source models because it is capable of bounding the output parameter range and will restart from a new initial simplex if a false valley is reached until a successful convergence is achieved. The optimization was performed using the full-waveform data filtered with a Butterworth bandpass filter between 0.2-7 Hz. Due to the variable duration of a chugging sequence, the portion of the signal that is used is adjusted to encompass the entire chugging signal while omitting the initial explosion. Only infrasound data were used in this analysis.
Nelder-Mead is a commonly used optimization technique because it is comparatively time-efficient relative to other methods such as the Spendley method or many global convergence optimizers. The bounded NM function nmkb from the R package dfoptim [Varadhan et al. 2018] is used with a RMS objective function calculated as the difference between the output model synthetic and recorded tremor data.
Although both the three-layered conduit model [Garcés 1997] and the oscillating gas pocket model [Girona et al. 2019] require several input parameters to calculate the theoretical result, certain parameters with characteristic values remain fixed during the simulation, while others are allowed to vary within the constraints permitted in the NM optimization. Due to the fundamental differences in structure between the two models, there is an imbalance of variable input parameters in which the layered-conduit model accepts six variables while the gas pocket model only accepts three. While this discrepancy is a possible source of bias in the comparison of the models, the spectral and waveform results can still be compared with the original chugging sequence to assess the viability of each model for chugging analysis.
Variable parameters fed to the NM algorithm in the layered-conduit model are the conduit section lengths L i and sound speeds c i to account for one geometric and one fluid compositional parameter. Other parameters such as magma density (ρ i ) and viscosity (µ i ) are fixed during the simulation with values typical of andesitic and basaltic magma compositions [Chesner and Rose 1984;Eichelberger and Izbekov 2000;Hanson et al. 2010;Webb and Dingwell 1990]. Variable parameters iterated during gas pocket model NM optimization are the thickness of the permeable cap (L), the permeability of the cap (κ), and the initial thickness of the gas pocket (D) prior to expansion due to gas impact.

Results
The model output parameters determined with NM optimization show that comparatively little variance is observed across output parameters related to the melt composition in the oscillating gas pocket model, whereas the standard deviation of the speed of sound parameter in each conduit layer either exceeds or approaches the boundaries initially set at the beginning of the NM algorithm (Tables 1 and 2). Parameters related to the geometry of the system also generally exhibit lower variance in the gas pocket model. Comparisons of the output parameters of the two models are difficult, however, due to differences in the model geometries and magma properties.
A summary of the constrained parameters calculated for 12 Karymsky records of chugging, 30 from Tungurahua, and 33 recorded at Fuego is provided in Tables 1 and 2. Fixed parameters used during the NM iterations for the layered conduit model include the fluid viscosity µ i " t1.2ˆ10´7; 1.4ˆ10´8; 1.6ˆ10´8u Pa s, magma density ρ i " t2450; 2500; 2700u kg m´3, and cross-sectional area of each layer S i " t10; 10; 20u m 2 . Non-variable parameters used in the gas pocket model include the gas viscosity µ g " 10´5 Pa s, temperature 1273.15 K, cap porosity φ " 10´4 %, and conduit radius R " 25 m.
Variable parameters fed to the NM algorithm are the conduit section lengths L i and sound speeds c i . Allowed ranges for the sound speed in the upper two conduit sections are considerably lower than the propagation velocity in the atmosphere, estimated to be about 320 m s´1. Two-phase magma flow regimes are highly sensitive to changes in gas abundance [McWilliam and Duggins 1969]. Phase changes from pure liquid to a two-phase flow can cause order of magnitude changes in sound velocity in magma; increase in isentropic com-pressibility as a function of 1{p due to the addition of gas, as well as the reduced density of the liquid magma, causes the velocity to decrease below the expected pure liquid and gas phase magma velocities.
Sound speeds in multiphase flow regimes begin to decrease once exceeding temperatures around 373 K in both pure and mixed phases, and will continue to reduce below characteristic sound speeds in pure gas phases, thus justifying the sound speed ranges (c 1,2 "10-250 m s´1, c 3 "1000-2500 m s´1) used in the NM optimization [Kieffer 1977]. The constrained variable parameters in the gas pocket model are the thickness of the permeable cap (L "2-30 m), the permeability of the cap (κ " 5ˆ10´9-2ˆ10´8 m 2 ), and the initial thickness of the gas pocket (D " 0.03-0.5 m). These initial ranges were the same for NM calculations performed on chugging data recorded at all three volcanoes to reduce variation in optimized parameter outputs.
Resolution of the spectral attributes of chugging data is a primary indication of the model's ability to estimate the source characteristics of a chugging sequence. Figure 5 shows the spectral output from both models for a single chugging sequence recorded at each volcano. A vertical red dashed line has been added to both modeled spectra to indicate the peak frequency of the chugging signal. The gas pocket model accurately predicts the 0.6-1.4 Hz dominant frequency observed in the data, as well as some secondary dominant frequencies and high frequency noise. The layered conduit model is in some cases able to resolve the dominant chugging frequency, but in others the predicted dominant frequency aligns with a secondary peak observed in the data. Some synthetic signals produced by the layered conduit model fail entirely to reproduce the amplitude variations in the time series data or the calculated dominant frequency does not align with any observed frequency peak. The modal structure in this model is extremely dominant and suppresses interpeak and high frequency fluctuations. Neither model resolves the characteristic asymmetry often observed in chugging signals, although the gas pocket model appears to possibly resolve double-periodicity observed in some sequences.

Discussion
The purpose of this research is to determine a set of reasonable source characteristics that effectively describe the spectral and physical properties of chugging sequences recorded at Karymksy, Tungurahua, and Fuego volcanoes, as well as to consider the validity of two general harmonic tremor models for understanding chugging mechanisms. A comparison of the two models used in this work reveals that the spectral characteristics of chugging are better resolved by a model that describes the mechanism driving chugging signals as a

Presses universitaires de �rasbourg
Page 257 Modeling source parameters of quasi-periodic tremor

Traphagan & Lees, 2020
Table 1 -Output parameters from the three-layered conduit model [Garcés 1997]. Parameters were calculated using an iterative Nelder-Mead algorithm. The values shown are averaged across 12 records from Karymsky, 30 at Tungurahua, and 33 at Fuego. shallow oscillating gas pocket trapped beneath a permeable cap, rather than a conical layered resonant conduit. Although dominant frequencies are often resolved by the resonant conduit model, in all model trials, high frequency fluctuations and frequency amplitudes observed in the data are not well constrained. However, the gas pocket model [Girona et al. 2019] does resolve the dominant frequencies and maintain a similar spectral structure in the high frequency bands to the recorded signal ( Figure 5). Amplitude envelope and double-periodicity are also resolved by this model, as shown in Figure 6. Thus, we suggest that the oscillating gas pocket model provides an adequate estimate of the source dynamics driving chugging. The physical interpretation from this model of the conduit suggests that gases may be initially released from a "16-20 m deep magma reservoir due to an explosion or earthquake, which corresponds to the initial high-frequency signal detected prior to a chugging episode shown in Figure 1. The exsolved gases may then travel through a conduit in packages or clouds until impacting and expanding the gas pocket. Grouped impacts provide an explanation for the emergent nature of chugging signals, while a series of multiple rising gas packages could generate intermittent tremor sequences as observed in some chugging signals recorded at Tungurahua and Karymksy.
The transition time between the gases being released from the reservoir and contacting the gas pocket also explains the time-lag observed in all chugging signals. Using an understanding of the pressure and fluid dynamics of the magma through which the gases travel before impact, a model could be developed in the future to estimate the conduit length between the magma reservoir and the bottom of the plug.
The physical representations of the gas pocket model are supported by geophysical observations at several gas releasing volcanoes. For example, chemical and temperature data from the crater lake at Ruapehu volcano, New Zealand, suggest the presence of a singlephase gas or vapor region exsolving from a lower magma solution [Christenson and Wood 1993]. Seismic analysis of harmonic tremor recorded at Ruapehu by Hurst and Sherburn [1993] suggests that the dominant 2 Hz frequency of the signal is explained by a resonator in the same gas region as proposed by Christenson and Wood [1993]. The concept of gas pocket expansion due to pressurization under a conduit plug is also proposed as a model for generating harmonic tremor following earthquake swarms or explosive eruptions at Sakurajima volcano, Japan [Maryanto et al. 2008], as well as for regular gas venting and tilt (inflation/deflation) cycles observed at Santiaguito volcano, Guatemala ]. These conduit structures have been Presses universitaires de �rasbourg  [Girona et al. 2019]. Results for other Tungurahua, Karymsky, or Fuego sequences can be provided on request. Normalized amplitude spectra are the same as the Tungurahua results in Figure 5. Data was recorded on May 28, 2010.
inferred to provide a physical coupling of shallow resonance supported by seismic data and continuous surface gas release observed at each volcano. Chugging signals have also been suggested to originate from a shallow subsurface oscillation of accumulated gases in Johnson and Lees [2000] and Lees et al. [2004].
The parameter ranges used to estimate chugging source dynamics are constrained by petrologic and seismic analysis of explosions leading to the formation of a conduit plug [Iverson et al. 2006]. Low surrounding country rock porosity (φ ď 15 %) and high magma flux (Q ě75 m 3 s´1) at z ď 100 m have been suggested as favorable conditions for the formation of a thin (L ď 60 m) solid permeable cap [Diller et al. 2006], consistent with the input ranges and estimated parameters modeled in this research. While the shallow gas pocket model [Girona et al. 2019] provides a sufficient method for estimating many characteristics of chugging and produces output parameters consistent with observations of typical plugged volcanoes, we do not discount the layered conduit model for harmonic tremor generated by deeper or larger conduit systems. Within the scope of this research, we suggest that the plugged conduit model is capable of resolving many recorded properties of quasi-periodic tremor.
Several alternative models have been proposed to describe the spectral and temporal characteristics ob-served from seismoacoustic harmonic tremor signals. Notably, Chouet [1985] models the source of tremor seismicity as a vertically buried cylindrical resonator coupled to a homogeneous halfspace. Julian [1994] models tremor-generating processes analogously to the resonance in a wind instrument, where nonlinear oscillations are excited into motion due to the flow of an incompressible Newtonian fluid through a deep viscoelastic channel. While these models are able to account for many of the characteristics observed in some harmonic tremor signals, the signal processing is strictly adjusted to analyze seismic signatures only. Lees and Bolton [1998] proposes similar methods to Julian [1994] for chugging-modeling and theoretically can resolve many chugging features, although the full development of this model remains a future objective. The approach in Lees and Bolton [1998] is conceptually similar to the mechanics of a venting pressure cooker, in which gas pressure is continuously balanced by an unstable equilibrium between the external downward force of the cap and the internal pressure of the body due to a constant flux of gas from a deeper conduit system.
Experimental methods have also been applied to describe harmonic tremor as an outgassing flux-driven process, such as Hellweg [2000], while other models have described the source of harmonic tremor Traphagan & Lees, 2020 as a horizontal swaying motion of a magma column [Bercovici et al. 2013;Jellinek and Bercovici 2011], repetitive pressure transients from bubble bursting [Lesage et al. 2006], and a repeating exploding source in a resonant conduit [Garcés and McNutt 1997]. In consideration of these models, we selected the models derived in Girona et al. [2019] and Garcés [1997] to provide a comparison of two infrasound-focused models that could be used to resolve the major characteristics of chugging.

Conclusions
This research adapted two harmonic tremor models from previous works to estimate the source properties of quasi-periodic harmonic tremor signals recorded at Karymksy, Tungurahua, and Fuego volcanoes. The model developed in Garcés [1997] describes the volcanic magma conduit conically, increasing in width with depth. The model is geometrically structured with three stacked layers, each with a unique set of impedance properties and dimensions, through which longitudinal standing waves are driven by the injection of fluid from a deeper magma body and propagated into the atmosphere via an open vent at the top. The second model used in this study was originally developed by Girona et al. [2019], which models the oscillation of pressure due to gas bubble bursting upon impact with a gas pocket trapped beneath a permeable overlying cap.
We suggest that the source mechanisms of chugging are best explained by the oscillating gas pocket model by Girona et al. [2019] due to the successful resolution of the dominant chugging frequency of 0.6-1.4Hz and low variance observed in the output parameter values. The model is also consistent with direct observations made at these volcanoes that chugging events typically occur simultaneously with degassing and shallow conduit processes. The model predicts an average cap thickness, initial gas pocket thickness, and cap permeability across all three volcanoes of approximately 19.6 m, 0.05 m, and 5.0ˆ10´9 m 2 . These fluid and conduit parameters provide an estimate of the physical dynamics at each volcano, but several processes not addressed by the oscillating gas pocket model could be occurring simultaneously to create the complex features observed in chugging signals.
Although the oscillating gas pocket model successfully recovered many chugging characteristics, there are several properties of chugging signals that remain unexplained. For example, the characteristic asymmetry of many signals observed in the acoustic data were not resolved by either model used in this study. Therefore, a combination of multiple models could be used to more accurately describe the conduit dynamics driving this process.