Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Dehghanniri, Vahid
Master of Science - MSc
University of British Columbia
Pre-eruptive seismic tremor with similar spectral properties is observed at active volcanoes with widely ranging conduit geometries and structures. Accordingly, the ``magma wagging'' model introduced by Jellinek & Bercovici[21] and extended by Bercovici et al.[6] hypothesizes an underlying mechanism that is only weakly-sensitive to volcano architecture: Within active volcanic conduits, the flow of gas through a permeable foamy annulus of gas bubbles excites and maintains an oscillation of a central magma column through a well-known Bernoulli effect. In this thesis, we carry out a critical experimental test of this underlying mechanism for excitation. We explore the response of analog columns with prescribed elastic and linear damping properties to forced annular airflows. From high-speed video measurements of linear and orbital displacements and time series of accelerometer measurements we characterize and understand the excitation, evolution, and steady-state oscillating behaviors of analog magma columns over a broad range of conditions. We identify three distinct classes of wagging: i. rotational modes which confirms predictions for whirling modes by Liao et al.[26]; as well as newly-identified ii. mixed-mode; and iii. chaotic modes. We find that rotational modes are favored for symmetric, and high intensity forcing. Mixed-mode responses are favored for a symmetric and intermediate intensity forcing. Chaotic modes occur in asymmetric or low intensity forcing. To confirm and better understand our laboratory results, and also extend them to conditions beyond what is possible in the laboratory, we carry out complementary two-dimensional simulations of our analog experiments. Our combined experimental and numerical results can be applied to make qualitative predictions for natural testable in future studies of pre- and syn-eruptive volcano seismicity. Long before an eruptive phase, the total gas flux is low and we expect magma wagging in a chaotic mode, independent of the spatial distribution of the gas flux. At a pre-eruptive state signaled by gas flux increasing, if the distribution of gas flux is approximately symmetric, we expect a transition to mixed and possibly rotational wagging modes. During an eruption, fragmentation and explosions within the foamy annulus can cause spatial heterogeneity in permeability resulting in non-uniform gas flux that favors chaotic wagging behavior.
An Experimental Study of Volcanic TremorDriven by Magma WaggingbyVahid DehghanniriBSc, University of Tehran, 2010MSc, K.N. Toosi University of Technology, 2013a thesis submitted in partial fulfillmentof the requirements for the degree ofMaster of Scienceinthe faculty of graduate and postdoctoralstudies(Geophysics)The University of British Columbia(Vancouver)January 2020c© Vahid Dehghanniri, 2020The following individuals certify that they have read, and recommend tothe Faculty of Graduate and Postdoctoral Studies for acceptance, the thesisentitled:An Experimental Study of Volcanic Tremor Driven by MagmaWaggingsubmitted by Vahid Dehghanniri in partial fulfillment of the requirementsfor the degree of Master of Science in Geophysics.Examining Committee:A. Mark JellinekSupervisorKelly RussellSupervisory Committee MemberValentina RadicSupervisory Committee MemberAli AmeliAdditional ExamineriiAbstractPre-eruptive seismic tremor with similar spectral properties is observed atactive volcanoes with widely ranging conduit geometries and structures. Ac-cordingly, the “magma wagging” model introduced by Jellinek & Bercovici[21]and extended by Bercovici et al.[6] hypothesizes an underlying mechanismthat is only weakly-sensitive to volcano architecture: Within active volcanicconduits, the flow of gas through a permeable foamy annulus of gas bubblesexcites and maintains an oscillation of a central magma column through awell-known Bernoulli effect. In this thesis, we carry out a critical exper-imental test of this underlying mechanism for excitation. We explore theresponse of analog columns with prescribed elastic and linear damping prop-erties to forced air annular airflows. From high speed video measurements oflinear and orbital displacements and time series of accelerometer measure-ments we characterize and understand the excitation, evolution, and steady-state oscillating behaviors of analog magma columns over a broad range ofconditions. We identify three distinct classes of wagging: i. rotationalmodes which confirms predictions for whirling modes by Liao et al.[26]; aswell as newly-identified ii. mixed-mode; and iii. chaotic modes. We findthat rotational modes are favored for symmetric, and high intensity forcing.Mixed-mode responses are favored for a symmetric and intermediate inten-sity forcing. Chaotic modes occur in asymmetric or low intensity forcing.To confirm and better understand our laboratory results, and also extendthem to conditions beyond what is possible in the laboratory, we carry outcomplementary two-dimensional simulations of our analog experiments.Our combined experimental and numerical results can be applied to makeiiiqualitative predictions for natural testable in future studies of pre- and syn-eruptive volcano seismicity. Long before an eruptive phase, the total gas fluxis low and we expect magma wagging in a chaotic mode, independent of thespatial distribution of the gas flux. At a pre-eruptive state signaled by gasflux increasing, if the distribution of gas flux is approximately symmetric,we expect a transition to mixed and possibly rotational wagging modes.During an eruption, fragmentation and explosions within the foamy annuluscan cause spatial heterogeneity in permeability resulting in non-uniform gasflux that favors chaotic wagging behavior.ivLay SummaryVolcanic tremor with a period of about 1 s is a common feature of mostexplosive eruptions and a major component of volcanic forecasting algo-rithms. “Magma wagging” model predicts that a column of magma canoscillate within a thin annular layer of foamy magma in a volcanic conduit,producing the observed seismic phenomena. Differential gas flow throughthis foamy annulus gives rise to pressure forces that drive and maintainthis oscillation continuously over hours to months. In this thesis, we useexperiments and numerical simulations to confirm the basic wagging phe-nomena. We identify new classes of magma column oscillation, dependingof the intensity and distribution of gas flows through the annulus. Our re-sults predict an evolution in the character and spatial correlation of tremorproperties that can be applied to eruption forecasting.vPrefaceChapter 1. Figures 1.2, 1.4, 1.5, 1.6, 1.7, 1.9, 1.10, 1.12, and 1.13 are usedwith permission from applicable sources.This thesis is ultimately based on the experimental apparatus and nu-merical modeling. The hardware design in Chapter 2 was done primarily byM. Jellinek and myself. The construction and tests were performed by S.Tomassi and I. The data analysis in Chapter 2 and Chapter 3 and numericalmodeling in Chapter 4 are my original work.A version of this material has been disseminated in form of poster presen-tation at 27th International Union of Geodesy and Geophysics Conference,Montreal, CanadaviTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . xxvDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xxvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Volcanic tremor . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Volcanic tremor source mechanism . . . . . . . . . . . . . . . 21.2.1 Previously proposed mechanisms . . . . . . . . . . . . 41.3 Magma wagging theory . . . . . . . . . . . . . . . . . . . . . 41.3.1 Magma wagging mechanism . . . . . . . . . . . . . . . 51.3.2 2D magma wagging with an impermeable annulus . . 71.3.3 2D magma wagging with a permeable annulus and gasflux variations . . . . . . . . . . . . . . . . . . . . . . 101.3.4 3D magma wagging model . . . . . . . . . . . . . . . . 17vii1.3.5 Summary of magma wagging model predictions . . . . 201.4 Critical knowledge gaps . . . . . . . . . . . . . . . . . . . . . 221.5 Thesis objectives and outlines . . . . . . . . . . . . . . . . . . 232 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 252.2 Column calibration . . . . . . . . . . . . . . . . . . . . . . . . 272.3 Calibrating the column response properties with mathemati-cal models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.1 2D mass-spring-damper model . . . . . . . . . . . . . 312.3.2 Cantilever beam model . . . . . . . . . . . . . . . . . 372.4 Scaling considerations . . . . . . . . . . . . . . . . . . . . . . 422.5 Energetic contributions to column top trajectories . . . . . . 453 Experimental Results and Discussion . . . . . . . . . . . . . 473.1 Rotational, mixed-mode, and chaotic column responses . . . . 473.2 Energy distribution . . . . . . . . . . . . . . . . . . . . . . . . 553.3 Spectral analysis . . . . . . . . . . . . . . . . . . . . . . . . . 593.4 Regime diagram . . . . . . . . . . . . . . . . . . . . . . . . . 603.5 Discussion of experimental results . . . . . . . . . . . . . . . . 653.5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . 663.5.2 Effects of forcing eccentricity e on column response . . 673.5.3 Effects of forcing intensity De on column response . . 693.5.4 Some unexpected results in e−De space . . . . . . . 734 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . 744.1 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . 754.1.1 Geometry and parameters . . . . . . . . . . . . . . . . 754.1.2 Equations and dimensionless parameters . . . . . . . . 754.1.3 Boundary conditions . . . . . . . . . . . . . . . . . . . 794.1.4 Mesh selection . . . . . . . . . . . . . . . . . . . . . . 804.1.5 Solvers and convergence . . . . . . . . . . . . . . . . . 804.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.2.1 A case study . . . . . . . . . . . . . . . . . . . . . . . 81viii4.2.2 Parametric study . . . . . . . . . . . . . . . . . . . . . 854.3 Discussions and summary . . . . . . . . . . . . . . . . . . . . 934.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . 934.3.2 The relation between the 2D numeric solution and ex-periments in terms of different class of responses: . . . 955 Implication for Real Volcanoes . . . . . . . . . . . . . . . . . 966 Concluding Remarks and Future Research Directions . . . 1006.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006.2 Limitations of the experiments . . . . . . . . . . . . . . . . . 1036.3 Future research directions: exploring volcano seismic data . . 104Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106ixList of TablesTable 2.1 Physical properties of VytaFlexTM columns used in theexperiments. ? comes from experiments, see §2.3.2 . . . . . 26Table 2.2 Values of damping coefficients and its standard deviationfor different columns. . . . . . . . . . . . . . . . . . . . . . 34Table 2.3 Values of fundamental damped frequency and its standarddeviation for different columns. . . . . . . . . . . . . . . . 34Table 2.4 Values of relaxation time, natural frequency of the system,and column’s characteristics number. . . . . . . . . . . . . 37Table 3.1 The average and standard deviation of time-averaged en-ergy components for all experiments in each response. . . . 57Table 3.2 The average value of energy ratio ER and its standarddeviation for different responses. . . . . . . . . . . . . . . . 59Table 4.1 List of parameters set in the model. . . . . . . . . . . . . . 77Table 4.2 List of non-dimensional parameters set corresponding toexperiments. . . . . . . . . . . . . . . . . . . . . . . . . . . 79xList of FiguresFigure 1.1 Time series of amplitude spectra prior to eruption on2009 −March − 28 03 : 24 at station REF, to the eastof Redoubt Volcano [20]. The spectrogram is plotted ona common logarithmic color scale using a 2-second 50%overlapping window. Time is aligned to the approximateexplosion onset at t = 0, and frequency is plotted upto the broadband Nyquist frequency of 25 Hz. Tremorduring transitions from background volcanic unrest toeruption begins as a harmonic or nearly monochromatic(‘narrow-band’) signal. As the eruption progress, boththe maximum frequency and the total signal bandwidthincrease to become ‘broadband’. Tremor gliding can beobserved in the last three minutes prior to eruption. . . . 3xiFigure 1.2 From [21] with permission. Sketch of the magma wag-ging model for volcanic tremor. Left panel: a very vis-cous magma column of radius Rm in a volcanic conduitof radius Rc is surrounded by a compressible permeablefoam annulus. The annulus is composed of stretched andcoalesced bubbles forming a tube-like matrix. Top rightpanel: an example of observed tube-like pumice [23]. Bot-tom right panel: an example of the observed gas ring atthe edge of Santiaguito volcano, which suggests the per-meable foam annulus. The depth H is typically of order1 km [8, 15]. . . . . . . . . . . . . . . . . . . . . . . . . . 6Figure 1.3 Sketch of 2D impermeable magma wagging model of vol-canic tremor showing an initial displacement to the right.The column is assumed to have a free surface at the topand is coupled viscously to the magma system below. As-sumed impermeable annulus prevents gas to rise relativeto annulus and plug of magma. Variables displaced in thediagram are discussed in the accompanying text. . . . . . 8Figure 1.4 From [21] with permission. Wagging frequency and shearstrain rate versus annulus thickness Rm/Rc for severalconduit radii Rc. Wagging frequency is calculated usingequation 1.2. Strain rate is (Q/(ρmpiR2c))/(Rc−Rm) andthe critical strain rate is (cG/µm), where Q is a typicalmass flow rate, G is the shear modulus of the melt, and cis a geometric constant. For more information see Fig. 3in [21]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10xiiFigure 1.5 From [6] with permission. Sketch of the cylindrical an-nulus for the permeable magma wagging model showinga displacement to the right. An important different toFig. 1.3 is that this displacement is driven as a result ofa pressure difference from left to right related to highergas speed on the right side (the Bernoulli’s effect). Thecolumn is assumed to have a free surface at the top and iscoupled viscously to the magma system below. Variablesdisplaced in the diagram are discussed in the accompany-ing text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 1.6 From [6] with permission. Oscillations in displacementu of the magma column at the column top, as predictedby the linear stability analysis. The oscillations are thesuperposition of only two modes, the fundamental modeand the least stable mode. Figures shown are for variousvalues of M, λ and η as indicated. All frames have thesame values of φ0 = 0.7, β = 100 and γ = 0. . . . . . . . 14Figure 1.7 From [6] with permission. Oscillations in displacement uof the magma column at the column top, as predicted bythe non-linear stability analysis. The oscillations are thesuperposition of only two modes, the fundamental modeand the least stable mode. Figures shown are for variousvalues of M, λ and η as indicated. Both frames have thesame values of β = 100 and γ = 2× 10−5. . . . . . . . . . 15Figure 1.8 Seismic and Acoustic time series for Karymsky and San-gay volcanoes [22] show growing sequence of beating en-velopes of tremor as well as correlated tremor and gas fluxvariations (acoustic time series). . . . . . . . . . . . . . . 16xiiiFigure 1.9 From [30] with permission. a. Time series for both 10-s averaged tremor amplitude, and SO2 vent gas emissionrate (in kg/s) for Fuego volcano in 2009. The top and bot-tom abscissa are offset by 32 s to compensate for overalllag. b. Evolution of the time lag between tremor and gasflux, where larger diamonds indicate correlation betweentime series larger than 0.65. . . . . . . . . . . . . . . . . . 17Figure 1.10 From [6] with permission. Time series and power spec-trum for surface displacement for a case with λ = 0.5,M = 0.1, and η = 1. The power spectrum is from aFourier analysis of the time series with the backgroundgrowth removed; power is in terms of U2(ω) where U(ω)is the discrete Fourier transform of U(H), and ω is thediscrete angular frequency. . . . . . . . . . . . . . . . . . 18Figure 1.11 Unbiased cross-correlation function between signals fromseismometers located north (RDN) and south (RSO) ofRedoubt volcano prior to its eruption on 2009 March 22.There is a positive correlation between the signals witha time lag of 0.55 s. The time lag is approximately halfthe period of the period of the dominant signal (1 s).Dash lines show the 95% confidence interval for cross-correlation calculations. . . . . . . . . . . . . . . . . . . . 19Figure 1.12 From [26] with permission. Sketch of the 3D magma wag-ging model. The column is assumed to have a free surfaceat the top and is coupled viscously to the magma systembelow. φ is the polar angle of the magma columns dis-placement, u is the radius displacement, or equivalentlythe displacement magnitude. Rc and Rm are the radii ofthe magma column and volcanic conduit, respectively. . . 20xivFigure 1.13 From [26] with permission. a. Contours of phase-shift ofP-waves generated by a 2D wagging at a frequency of 3.3Hz. Black lines and dashed lines are two groups of equal-phase contours which are out-of-phase with each other.Pentagons indicate the location of five virtual seismic sta-tions. The orange dash line with arrowhead indicates thedirection of the wagging plane.b. Contours of phase-shiftof P-waves generated by a counter-clockwise wagging withcircular whirling orbit at a frequency of 3.3 Hz. Blacksolid lines and blue dashed lines are two groups of equal-phase contours which are out-of-phase with each other.Pentagons indicate the location of five virtual seismic sta-tions. The orange dash line with arrowhead indicates thecounter-clockwise wagging direction. c. P-waves calcu-lated at the five virtual stations in a, which are in-phase orout-of-phase with each other. Stations that are not in thesame phase have half a period of oscillation T0/2 = 0.15 slag relative to each other.d. P-waves calculated at thefive virtual stations in b. . . . . . . . . . . . . . . . . . . 21Figure 2.1 Experimental setup and the arrangement of inlets in themanifold. Hollow cylinder dimensions: 10.16 cm internaldiameter × 180 cm long. Surrounding box dimensions:20 cm×20 cm×125 cm. The top camera is fixed 20 cmabove the apparatus, and the distance between side cam-eras and the tank surface is 150 cm. Distance betweenthe center of the manifold and inlets r is 4 cm, and inletsdiameter d is 0.95 cm. . . . . . . . . . . . . . . . . . . . . 28xvFigure 2.2 a. Example of trajectory of top of the column. Displace-ments in x and y directions are scaled by δ. Colorbarshows the dimensionless time t/Td. b. The correspond-ing time series in x direction, and c. The correspondingtime series in y direction. Transient-startup, steady-state,and damping regimes are defined in the text. . . . . . . . 29Figure 2.3 a. Example of trajectory of top of the column duringthe transient forcing, b. the corresponding time seriesin x direction, and c. the corresponding time series in ydirection. d. Example of trajectory of top of the columnduring the damping regime, e. the corresponding timeseries in x direction, and f. the corresponding time seriesin y direction. . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.4 Schematic of simplified version of experiments viewed fromtop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.5 a. The trajectory of top of the column in x direction, b.The trajectory of top the column in y direction. Dashedlines present envelopes fitted to trajectories to calculatedamping coefficient. . . . . . . . . . . . . . . . . . . . . . 34Figure 2.6 Histograms of damping coefficient (ζωn in 2.3) calculatedfrom damped time series. Solid lines represents the fittedGaussian curves. a., b., c., and d. correspond to C1, C2,C3, and C4 respectively. . . . . . . . . . . . . . . . . . . . 35Figure 2.7 a. power spectrum of Fig. 2.5 a., b. power spectrum ofFig. 2.5 b. . . . . . . . . . . . . . . . . . . . . . . . . . . . 36Figure 2.8 Histograms of damped angular frequency (ωd) calculatedfrom damped time series. Solid lines represent the fittedGaussian curves. a., b., c., and d. correspond to C1, C2,C3, and C4 respectively. . . . . . . . . . . . . . . . . . . . 36Figure 2.9 Schematic of simplified version of experiments from sideview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37xviFigure 2.10 Free body diagram for an element of column. w is exter-nally distributed load, V is shear force, and M is bendingmoment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 2.11 Different elliptical whirling orbits corresponding to differ-ent ER = KELin/KERot values. Circular orbit corre-sponds to ER = 0. When ER → ∞ orbit changes to astraight line. . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 3.1 Rotational response; a. The trajectory of top of thecolumn. Displacements in x and y directions are scaledby δ. The colorbar is the dimensionless time t/Td. b.Snapshots of column deflection during an orbit. c. Thecorresponding time series in x direction, and d. The cor-responding time series in y direction. Blue lines are upperand lower peak envelops of displacement time series. . . 49Figure 3.2 Trajectory of top of the column in a rotational response.Here flow is evenly distributed around the column e = 0and Deborah number is De = 1.4 × 104. Each subplotshows the trajectory during one dimensionless time Td,and the colorbar shows the corresponding time. The col-umn follows approximately constant orbits 2 to 3 times ineach time window. Trajectories are similar to each otherin each subplot, and the ones in other subplots in termsof radius and aspect ratio. The direction of rotation isconstant, i.e. it moves clockwise. . . . . . . . . . . . . . . 50Figure 3.3 Mixed-mode response; a. The trajectory of top of thecolumn. Displacements in x and y directions are scaledby δ. The colorbar is the dimensionless time t/Td. b.Snapshots of column deflection during an orbit. c. Thecorresponding time series in x direction, and d. The cor-responding time series in y direction. Blue lines are upperand lower peak envelops of displacement time series. . . . 51xviiFigure 3.4 Trajectory of top of the column in a mixed-mode response.Flow distribution around the column is approximatelyeven e = 0.08 and Deborah number is De = 6.8 × 103.Each subplot shows the trajectory during one dimension-less time Td and colorbar shows the corresponding time.The column mainly follows elliptical and circular orbits,e.g. t = 25, 26, 33 are elliptical orbits while t = 34 and 36are circular. Aspect ratios and radius of ellipses and cir-cles are different from time to time, e.g. the evolution oftrajectory from t = 30 to t = 33 as the orbits are enlarg-ing. Sometimes we observe linear movements, e.g. t = 28.Although trajectories change over time but the directionof motion is always the same, i.e. for this example columnmoves in clockwise direction. . . . . . . . . . . . . . . . . 52Figure 3.5 Chaotic response; a. The trajectory of top of the col-umn. Displacements in x and y directions are scaled by δ.The colorbar is the dimensionless time t/Td. b. Snapshotsof column deflection during a cycle. c. The correspondingtime series in x direction, and d. The corresponding timeseries in y direction. Blue lines are upper and lower peakenvelops of displacement time series. . . . . . . . . . . . . 53xviiiFigure 3.6 Trajectory of top of the column in a chaotic response.Flow distribution around the column is highly unevene = 1 and Deborah number is De = 75. Each subplotshows the trajectory during one dimensionless time Tdand colorbar shows the corresponding time. Linear tra-jectories are the most observed pattern in this figure whilethe direction of linear trajectories are different, e.g. tra-jectory at t = 41 is parallel to x axis while at t = 46 isparallel to y axis. At t = 40 and t = 44 we see that di-rection of linear motion changes. There are two ellipticaltrajectories at t = 37 and t = 39 while the orientationand size of ellipses are different. Other trajectories area combination of linear and circular movements in smallscales as we see at t = 36 or t = 43. . . . . . . . . . . . . 54Figure 3.7 a. Energy time series in a rotational response. Flowconditions are e = 0 and De = 1.4 × 104., b. Energytime series in a mixed-mode response. Flow conditionsare e = 0.08 and De = 6.8 × 103., c. Energy time se-ries in a chaotic response. Flow conditions are e = 1 andDe = 75. Time is scaled with Td and rotational kineticenergy KERot, linear kinetic energy KELin, potential en-ergy PE, and total specific energy E0 are normalized bythe maximum of specific energy E0max during each exper-iment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Figure 3.8 Distribution of time-averaged specific energy componentsfor all the experiments in a. rotational responses, b.mixed-mode responses, and c. chaotic response. d. Dis-tribution of energy ratio ER for different motions. Dashedlines are fitted Gaussian curves to ER distributions, detailin Table 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . 58xixFigure 3.9 Example of Power spectrum for a. rotational, b. mixed-mode, and c. chaotic response. Red lines are from spec-tral analysis of accelerometers data and black lines comefrom spectral analysis of trajectory time series. The fre-quency is normalized by the damped frequency of the col-umn fd. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 3.10 Power spectrum of all the experiments in a. rotational, b.mixed-mode, and c. chaotic response. The solid black linerepresents the average value power spectrum for all theexperiments in each mode. The frequency is normalizedby the damped frequency of the column fd. . . . . . . . . 62Figure 3.11 Effects of forcing eccentricity e and intensity De on col-umn responses. Shapes represents responses as: © isrotational, ♦ is mixed-mode, and  is chaotic response.The colorbar shows the value of ER. . . . . . . . . . . . . 63Figure 3.12 Effects of forcing eccentricity e and intensity De on C1 re-sponses. Shapes represents responses as: © is rotational,♦ is mixed-mode, and  is chaotic response. The colorbarshows the value of ER. . . . . . . . . . . . . . . . . . . . 64Figure 3.13 Effects of forcing eccentricity e and intensity De on C2responses. Shapes represents responoses as: © is rota-tional, ♦ is mixed-mode, and  is chaotic response. Thecolorbar shows the value of ER. . . . . . . . . . . . . . . 64Figure 3.14 Effects of forcing eccentricity e and intensity De on C3 re-sponses. Shapes represents responses as: © is rotational,♦ is mixed-mode, and  is chaotic response. The colorbarshows the value of ER. . . . . . . . . . . . . . . . . . . . 65Figure 3.15 Effects of forcing eccentricity e and intensity De on C4 re-sponses. Shapes represents responses as: © is rotational,♦ is mixed-mode, and  is chaotic response. The colorbarshows the value of ER. . . . . . . . . . . . . . . . . . . . 66xxFigure 3.16 Schematic of gas velocity and pressure variation aroundrubber column, and imparted pressure at the inlet for asymmetric airflow condition. Longer arrows and darkercolors represent bigger magnitude of flow velocity andpressure forces. Note, imposed gas pressure force at baseacting on a deflected column is always symmetric. . . . . 70Figure 3.17 Schematic of gas velocity and pressure variation aroundrubber column, and imparted pressure at the inlet for anasymmetric airflow condition. Longer arrows and darkercolors represent bigger magnitude of flow velocity andpressure forces. Note, pressure variations imparted atbase are not symmetric. Even if the column would prefera rotational mode, this inlet condition always perturbsthe column away from axi-symmetry. . . . . . . . . . . . . 71Figure 4.1 Schematic of the model used in COMSOL Multiphysics.An elastic column in dimensions of D × L is fixed at thebottom. The initial gap width between the elastic columnand the surrounding wall is δ. Different types of meshused in this study, with a number of elements, extra finemesh: 93940, finer mesh: 22550, fine: 15654, and normal:9968. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 4.2 Top displacement of the column with 4 different meshsettings. In this model KG = 10, Re = 1000, m = 850,and χ = 0.5. The time required for each simulation was45, 125, 224, and 410 min for Normal, Fine, Finer, andExtra F ine mesh settings, respectively. . . . . . . . . . . 80xxiFigure 4.3 Column displacement at L/3, L/2, 2L/3, and L distancesfrom the fixed side where L is length of the elastic column.Arrows and dash lines show the period of low- and high-frequency oscillations. Non-dimensional parameters forthis simulation are KG = 0.1, Re = 100, m = 850, andχ = 0.5. y axis is scaled by gap width δ and time is scaledby T = L/V . . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 4.4 Average velocity of the airflow in right side and left side ofthe elastic column at a. L/3, b. 2L/3, and c. L distancesmeasured from the bottom of the channel. Arrows anddash lines show the period of low and high-frequency os-cillations. Parameters are KG = 0.1, Re = 100, m = 850,and χ = 0.5. Time is scaled by T = L/V . . . . . . . . . . 83Figure 4.5 Pressure variations of the airflow in right and left sides ofthe elastic column at a. L/3, b. 2L/3, and c. L distancesmeasured from the bottom of the channel. The solid anddashed black lines present the pressure envelopes at theright and left the side of the channel, respectively. Arrowsand straight dash lines show the period of low and high-frequency oscillations. . . . . . . . . . . . . . . . . . . . . 84Figure 4.6 a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the rightchannel at 2L/3 from the bottom of the channel. Dashedlines in b. represent velocity at the left channel. Arrowsand dash lines show the period of low and high-frequencyoscillations. Non-dimensional parameters for these simu-lations are KG = 1, Re = 1000, m = 850, and χ changesfor each simulation. Time is scaled by T = L/V . . . . . . 86xxiiFigure 4.7 a.Top displacement of the column, b. average flow veloc-ity, and c. pressure variation of the flow in the right of thechannel at 2L/3 from the bottom of the channel. Parame-ters are KG = 1, χ = 0.1, m = 1000, and for different Re.Arrows and dash lines show the period of low and high-frequency oscillations. Flow velocity is magnified 100×and 10× for Re = 10 and Re = 100, respectively. Pres-sure variations for Re = 10 and Re = 100 is 104× and100× magnified, respectively. . . . . . . . . . . . . . . . . 88Figure 4.8 a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the rightchannel at 2L/3 from the bottom of the channel. Arrowsand dash lines show the period of low and high-frequencyoscillations. Non-dimensional parameters for these sim-ulations are χ = 0.1, Re = 1000, m = 1000, and KGchanges for each simulation. Time is scaled by T = L/V . 90Figure 4.9 a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the rightchannel at 2L/3 from the bottom of the channel. Non-dimensional parameters for these simulations are χ = 0.9,Re = 1, m = 100, and KG = 0.01. Time is scaled byT = L/V . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 4.10 a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the rightchannel at 2L/3 from the bottom of the channel. Arrowsand dash lines show the period of low and high-frequencyoscillations. Non-dimensional parameters for these simu-lations are KG = 1, Re = 1000, χ = 0.5, and m changesfor each simulation. Time is scaled by T = L/V . . . . . . 92Figure 4.11 Top displacement of the column. Non-dimensional pa-rameters for this simulation are KG = 0.01, Re = 1,χ = 0.9, and m = 100. . . . . . . . . . . . . . . . . . . . . 95xxiiiFigure 5.1 Some snapshots of observed gas ring at the edge of Karym-sky volcano. . . . . . . . . . . . . . . . . . . . . . . . . . . 97Figure 5.2 Volcano’s activity in eccentricity e, Deborah De numberregime. A long time before the eruption, (t texplosion,volcanic degassing is inefficient, De ≈ 0, and degassingis uniform, e ≈ 0. Closer to eruption time, degassingenhances, De ↑, while annulus structure does not change.Eventually, enhancing degassing causes bubbles strain rateto increases above a critical strain rate, and fragmenta-tion and explosion occur. Fragmentation and explosioncause annulus to become heterogeneous and lead to non-uniform degassing around the volcano, e → 1. Each col-ored area corresponds to the ideal condition for one of therotational, mixed-mode, or chaotic regimes based on theresults in Chapter 3. . . . . . . . . . . . . . . . . . . . . 99xxivAcknowledgmentsFirst and foremost, I would like to express my sincere gratitude and ap-preciation to my mentor and supervisor Prof. Mark Jellinek. Mark thanksfor being so patient toward my painfully-naive questions. Thanks for guid-ing me to become a scientist and importantly being so knowledgeable andcaring.I would also like to thank the members of my advisory committee: Prof.Kelly Russell and Prof. Valentina Radic. Kelly! thanks for introducing thebeautiful world of volcanoes to me. Valentina! thanks for showing me howto create something meaningful out of big bunch of data.I would also like to thank Jorn Unger and Santiago Tomassi. Everytime I had a problem with the experimental setup I came to Jorn and hehelped me patiently. Santiago, Thanks for helping me to run more than1000 experiments.I would also like to thank my wonderful lab-mates in MJCJ researchgroup for making my graduate life more enjoyable and fun.I am especially grateful to my parents. I thank them for their endlesssupport and unconditional love. This journey could have not been possiblewithout their support. I also thank my siblings for encouraging me to followmy passion. I am thankful to all of them and forever grateful for havingthem in my life.I would also like to thank my beautiful wife, Parisa. If it was not foryou, I would not stand at this point.Last but certainly not least, Id like to thank my new family for support-ing me every step of the way.xxvDedicationTo my parentsfor their unconditional loveChapter 1IntroductionViolent volcanic eruptions such as the 1980− 86 Mount Saint Helens event[12], the 1991 eruption of Mount Pinatubo [19], the 1995− 2013 eruption ofSoufriere Hills [38], or the 2010 Eyjafjallajkull [17] event can be devastatingto local populations and can disrupt air travel with profound economic con-sequences felt worldwide. With the current steady rise in the populationsof cities on the flanks of active volcanoes and a growing reliance on interna-tional air travel to support an ever-increasing global economy, the need toreliably forecast catastrophic eruptions is paramount.The rest of this chapter is organized as follows. We initially introducevolcanic tremors as a promising indicator of a volcanic eruption. Next, webriefly discuss existing models of the source mechanism of volcanic tremorsand explain the reason that a universal model is needed. Later we explorethe recently proposed “magma wagging” model and discuss main phenom-ena that magma wagging predicts carefully. Then, we address some of thecavities in the literature of magma wagging and knowledge gaps that areneeded to be studied, and finally, we present our method to study magmawagging theory and improve it.11.1 Volcanic tremorAn important tool for establishing dangerous impending volcanic unrest ismonitoring of the volcanos seismic activity. A particular type of seismic ac-tivity called volcanic tremor almost always precedes the explosive volcaniceruptions that present the greatest danger to humans. Tremors persists forminutes to weeks and its seismic waveforms are characterized by a remark-ably narrow band of frequencies from about 0.5 Hz to 7 Hz [24, 28, 29, 31].Before major eruptions, tremor can occur in concert with increased gas fluxand related ground deformation [12, 27, 35]. Tremor frequency can varyin time. An intriguing feature of tremor observed at several eruptions is“gliding” within the frequencies of evenly spaced spectral peaks vary sys-tematically with time [20, 28] (see Fig. 1.1). McNutt and Nishimura [29]also show that around 90% of volcanic tremor decay exponentially at theend of eruptions. Also, some studies use volcanic tremor to estimate the Vol-canic Explosivity Index (VEI), a numeric scale that measures the relativeexplosivity of eruptions, of eruption while it is in progress [1]. Consequently,volcanic tremor forms a major component of eruption forecasting algorithms[7, 35], and understanding its origin might be a great help to forecast explo-sive eruptions.1.2 Volcanic tremor source mechanismAlthough tremor is an empirically useful indicator of imminent volcanism,a clear understanding of the underlying physical process is elusive. Severalstudies indicate that volcanic tremors have a common excitation mechanismwith a source that is repeatable, non-destructive, and at a fixed location [31].Benoit & McNutt [5] have shown evidence of changing apparent source lo-cation in space and suggest that tremor source is linear in shape rather thana point source. Unglert and Jellinek [37] identify three main classes of vol-canic tremor and discuss that there should be a generic relationship betweenspectral character and magma volcano characteristics, including magma vis-cosity, magma storage depth, and the physical properties of volcano edifices.2𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺𝐺`Redoubt station locationsFigure 1.1: Time series of amplitude spectra prior to eruption on 2009 −March − 28 03 : 24 at stationREF, to the east of Redoubt Volcano [20]. The spectrogram is plotted on a common logarithmic colorscale using a 2-second 50% overlapping window. Time is aligned to the approximate explosion onsetat t = 0, and frequency is plotted up to the broadband Nyquist frequency of 25 Hz. Tremor duringtransitions from background volcanic unrest to eruption begins as a harmonic or nearly monochromatic(‘narrow-band’) signal. As the eruption progress, both the maximum frequency and the total signalbandwidth increase to become ‘broadband’. Tremor gliding can be observed in the last three minutesprior to eruption.31.2.1 Previously proposed mechanismsSeveral models exist to explain aspects of the spectral properties volcanictremor source mechanisms. For example, tremor can be a continuous oscil-lation that reflects either a natural mode of vibration excited in the volcanicconduit, analogous to an organ pipes response to airflow [10, 16]. Alter-natively, tremors are caused by the elastic response of existing cracks tomagma flow, i.e., the rapid opening and closing of fractures or “crack flap-ping” [7, 34]. Hellweg et al. [18] also propose an “eddy shedding” modelin which conduit wall reaction to pressure fluctuations caused by eddies,which form behind obstacles, lead to volcanic tremor. “Pressure-cooker”and “soda-bottle” models [18, 22] explain periodic degassing and periodicash release controlled by some valve mechanism might be a tremor cause.Other proposed tremor mechanisms include magma and water interactions,bubble growth and collapse in magma or water [25], and resonance of magmabody under the ground [9]. All of these models explain the monochromaticproperties of tremor frequency.These approaches to understanding tremor focus only on the oscillationand do not address time-dependent properties of frequency or amplitude oftremor. Also, they involve very different physical phenomena. Furthermore,their success relies on an unrealistic assumption that all volcanoes world-wide are similar in terms of the geometry, structure, and constitution ofvolcanic conduit walls, as well as the gas content of the erupting magma.In addition to the initial mechanical structure of volcanic conduits beingdistinct, the evolutions of the magma-conduit system during eruptions willbe different from one volcano to the next. Ultimately, it is surprising thattremor characteristics are so consistent among volcanoes worldwide, and theuniversality of tremor properties particularly related to explosive volcanismremains, thus, a major enigma[21].1.3 Magma wagging theoryIn a departure from the conventional focus on the geometric and structuralproperties of volcanic conduits to explain the physical process behind tremor,4recent studies show that volcanic tremor may be a natural consequence ofthe flow dynamics of an ascending magma. The so-called “magma wagging”theory is proposed in [21], which employs the contemporary view of silicicvolcanic systems that ascending viscous magmas with high concentrations ofcrystals and bubbles form a stiff columnar plug surrounded by a compressiblelayer of sheared foam (Fig. 1.2).1.3.1 Magma wagging mechanismIn silicic eruptions, involving high viscosity magma, ascent speeds of gasbubbles are negligible compared with the rise speed of magma [35]. As aresult, magma rises as a gas- and crystal-rich plug before and sometimesduring an eruption. Shearing at the conduit causes bubbles, which have amuch lower viscosity than the surrounding melt phase, to stretch, merge,and become interconnected [6]. Consequently, a vesicular, permeable, andcompressible annulus of bubbly magma forms around the magma plug. Mod-els [13, 32, 33, 35] as well as some observations of pumice with porosity φ0of 30 − 90% from effusive and explosive eruptions represents of permeabletube-like matrix within the annulus [13, 35], see Fig. 1.2.For most geologically-relevant conditions, the stiff magma column that isconnected to a magma reservoir below and has an upper free surface will os-cillate or “wag” laterally or azimuthally against a restoring gas-spring forceof the compressible foam layer. When the magma column is displaced fromits central resting position, the annulus is compressed in the direction ofcolumn displacement and is decompressed in the other direction. Throughthe ideal gas law, there will be net pressure force originated from permeablefoamy-springy annulus pushing back the column toward its resting position,provided that the column displacement occurs over a time scale that is smallcompared to the time for compressed gas to drain by porous media flow. Thecolumns inertia causes it to overshoot the resting position during an oscil-lation. Following an initial displacement, the column oscillates inside theannulus until viscous forces in the magma column damp it. Pressure forcesimparted at conduit walls result in seismic energy being transmitted through5Figure 1.2: From [21] with permission. Sketch of the magma waggingmodel for volcanic tremor. Left panel: a very viscous magmacolumn of radius Rm in a volcanic conduit of radius Rc is sur-rounded by a compressible permeable foam annulus. The an-nulus is composed of stretched and coalesced bubbles forminga tube-like matrix. Top right panel: an example of observedtube-like pumice [23]. Bottom right panel: an example of theobserved gas ring at the edge of Santiaguito volcano, which sug-gests the permeable foam annulus. The depth H is typically oforder 1 km [8, 15].6the conduit walls surrounding the wagging magma column and leading tothe tremor signal observed on local seismic stations [6, 21].1.3.2 2D magma wagging with an impermeable annulusIn the original model, Jellinek & Bercovici [21] make five simplifications.First, they assume that bubble- and crystal-rich magma, and the surround-ing annulus foam layer rise together. Second, they assume that inner columnof magma of radius Rm is cylindrically axisymmetric and initially at rest andcentered inside a cylindrical conduit of radius Rc. Thus, the resting gap be-tween the conduit wall and the inner column is L = Rc − Rm. Third, tofacilitate analysis, they define the (linear) column instability to lateral dis-placements in polar coordinates such that when a column is displaced leftor right by u(z), the maximum gap width is L + u and the minimum gapwidth is L − u (see Fig. 1.3). Fourth, the annulus is assumed to be imper-meable. Fifth, they neglect the elastic behavior of the column itself as themaxwell relaxation time is 10−1 s, for rhyolitic magmas with viscositiesof 105−109 Pa and shear moduli in the range of 1010−1011 Pa [2], is shorterthan typical periods of volcanic tremor.With these assumptions, the original equation of motion governing wag-ging is akin to that for a linearly damped harmonic oscillator [21]:∂2~u∂t2= −ω20~u+ νm∂3~u∂z2∂t, (1.1)where ~u is the displacement from column resting position, νm is the magmakinematic viscosity, z is the vertical direction along with the column height.The natural angular frequency [rad/s] of the wagging magma column ω0depends on a balance between gas pressure-force (the spring force) and thecolumns inertia and is:ω0 =√2ρ0C2φ0ρm(R2c −R2m), (1.2)where the ordinary frequency [Hz] is ω0/2pi. Here C is the speed of sound7𝑧𝑧𝑥𝑥Figure 1.3: Sketch of 2D impermeable magma wagging model of vol-canic tremor showing an initial displacement to the right. Thecolumn is assumed to have a free surface at the top and is cou-pled viscously to the magma system below. Assumed imperme-able annulus prevents gas to rise relative to annulus and plugof magma. Variables displaced in the diagram are discussed inthe accompanying text.in gas (water vapor at magmatic temperature), ρ0 and φ0 are the densityand volume fraction of the undisturbed gas in the annulus, ρm is the magmadensity, and Rc and Rm are the conduit and magma column radii, respec-tively. The left-hand side of (1.1) represents column acceleration. While onthe right-hand side, the first term is the acceleration related to the restoringfoam spring force in the annulus, and the second term is the decelerationarising in response to the production of vertical velocity gradient in themagma column as a result of viscous bending and wagging.The fundamental frequency of wagging (1.2) implies that increasing bub-ble volume fraction φ0 in the annulus decreases wagging frequency (ω0 ∝√1/φ0). Also, decreasing the column radius Rm for fixed Rc causes thewagging frequency to decrease (ω0 ∝√1/(R2c −R2m)). Critically and in8marked contrast to all other tremor models, ω0 is only very weakly sensitiveto conduit geometry through equation 1.2. Furthermore, this behavior leadsto an important prediction consistent with Fig. 1.1. InSAR monitoring [27]as well as seismic data [11] show that volcanos gas-pressurize prior to erup-tion. The cause of volcano pressurization is increased gas content inside themagma. A resulting increase in the bubble concentration in the foamy an-nulus reduces magma viscosity locally and enhances the localization of sheardeformation near the conduit walls [15], which leads to a narrower annulus,in turn. This progressive narrowing of the annulus width, L = Rc − Rm,causes the wagging frequency to climb.All the variables in equation 1.2 are well constrained except the conduitgeometry. Jellinek & Bercovici [21] show that using C = 700 m/s for a gaswhich is primarily water at about 1000 K, ρm/ρ0 = 100, 30% < φ0 < 90%,10 m≤ Rc ≤ 100 m, and 0.5 ≤ Rm/Rc < 1 the frequency of oscillationranges from 0.1 Hz to 5 Hz which is in the range of known volcanic tremors.Figure 1.4 shows wagging ordinary frequency ω0/2pi [Hz] on the left y-axis(solid lines) and shear strain rate in the annulus on the right y-axis (dashedlines) versus annulus thickness on x-axis for several conduit radii. Widercolumns wag in lower frequencies because of greater inertia forces they have.Figure 1.4 also shows the evolution in shear strain rate as annulus thick-ness L decrease with increased gas concentration. The shear strain rate can-not increase to infinity: Eventually, the shear strain rate reaches a maximumat which the time scale for the accumulation of shear strain becomes shorterthan time scale for viscoelastic relaxation of bubble walls, and sheared bub-bles rupture catastrophically and fragment [15, 36]. Thus, fragmentationis associated with a minimum annulus thickness and maximum frequency.The solid black line in Fig. 1.4 represents the threshold for fragmentation,and arrows correspond to minimum annulus thickness and maximum tremorfrequency.Fragmentation is envisioned to cause the annulus to evolve rapidly froma foam composed of tube-like bubbles with a homogeneous, on average, per-meability to a mixture of fractured tubes and shattered glass fragments witha permeability that is spatially heterogeneous [21]. A strong heterogeneity9Figure 1.4: From [21] with permission. Wagging frequency and shearstrain rate versus annulus thickness Rm/Rc for several conduitradii Rc. Wagging frequency is calculated using equation 1.2.Strain rate is (Q/(ρmpiR2c))/(Rc − Rm) and the critical strainrate is (cG/µm), where Q is a typical mass flow rate, G is theshear modulus of the melt, and c is a geometric constant. Formore information see Fig. 3 in [21].into the structure of the annulus causes gas volume fraction φ0 in the an-nulus to vary in space [21], causing a broadening of the wagging frequencybandwidth because ω0 ∝√1/φ0.1.3.3 2D magma wagging with a permeable annulus andgas flux variationsBecause the viscous resistance of the magma column to bending, (νm∂3~u∂z2∂t)in equation 1.1, will attenuate the oscillation over timescales as short ashours [6, 21], a driving mechanism to maintain the oscillation is necessary.Bercovici et al. [6] develop a hypothesis that the magma wagging model isa continuously self-excited system. They extend the magma wagging model10to include a permeable annulus, which allows gas to move relative to magmaand bubbles in the annulus (see Fig. 1.5) and allows for left-right variationsin gas flux.When the column is displaced, in the direction of displacement gas pres-sure increases and imposes a pressure force to restore the column displace-ment. At the same time, rise speed of compressed gas in the annulus in-creases through a nozzle effect. This increase in velocity on the compressedside causes a pressure drop which draws the column further in the direc-tion of displacement through the well-known Bernoulli effect. In contrast,gas decompresses in the opposite direction, its pressure decreases and leadto a pressure force to oppose the column displacement. At the same time,rise speed of decompressed gas in the annulus decreases. This decrease invelocity on the decompressed side cause a pressure increase which drawsthe column further in the direction of displacement. Therefore, differentialgas ascent speeds in annulus will induce lateral pressure variations that canexcite relatively low frequency “Bernoulli mode” that potentially suppliesenergy to higher frequency wagging.In more detail, as with the original model, Bercovici et al. [6] restrictwagging to be in 2D, in the x− z plane in Fig. 1.3. Thus, when the columnis displaced by u(z) the maximum gap width is L + u, and the minimumgap width is L − u. They consider gas flow in the annulus is driven bothby gas injected from below and modulated by left-right pressure variationsarising from lateral displacement of the magma column. They also considerthat the gas density, annulus porosity, and gas velocity are independent ofradial distance r from the center of the conduit.The greatest physical insight into 2D magma wagging model with perme-able annulus is via the the equations of mass and momentum conservationin dimensionless form:∂ρϕ∂t+∂ρϕw∂z= 0, (1.3)ρφ(∂w∂t+ w∂w∂z) = −∂ρ∂z− γ(β(1− φ) + ρφ), (1.4)11𝑧𝑧𝑥𝑥Figure 1.5: From [6] with permission. Sketch of the cylindrical an-nulus for the permeable magma wagging model showing a dis-placement to the right. An important different to Fig. 1.3 isthat this displacement is driven as a result of a pressure differ-ence from left to right related to higher gas speed on the rightside (the Bernoulli’s effect). The column is assumed to have afree surface at the top and is coupled viscously to the magmasystem below. Variables displaced in the diagram are discussedin the accompanying text.∂2u∂t2= − 2λβpi∫ 2pi0ρ(cosθ + 2λu cos(2θ))dθ + η∂3u∂z2∂t, (1.5)whereϕ = φ0 − u cosθ and φ = φ0 − (1− φ0)u cosθ,and the dimensionless control parameters are:γ =gUC2, β =ρmρo, η =νmCUand λ =U2Rm, (1.6)where U = (R2c − R2m)/2Rm and U ≈ L in the limit of L Rc. Hence,12γ represents the ratio of hydrostatic and dynamic gas pressure, β is theratio of magma and gas densities, λ expresses the annulus thickness to themagma column width, and the most important, η represents the ratio ofviscous forces arising in the magma column, depending on the rate of dis-placement, to the gas spring force in the annulus. Alternatively, η is theratio of time scale to build up elastic spring forces U/C to the time scale ofdamp oscillations through viscous effects U2/νm.To understand the response of magma column to lateral pressure varia-tions arising in response to gas flux variations, Bercovici et al. [6] investi-gate the linear and nonlinear stabilities of the 2D Cartesian magma waggingmodel (equations 1.3-1.5) analytically and numerically, respectively. At lin-ear order displacement u, gas velocity w, and density ρ all scale as eikz+st,where k is wavelength and s is the growth rate of a small perturbation. Theyfind the characteristics equation:s2 +φ0α2(s+ ikM)2k2 + φ0(s+ ikM)2+ ηk2s = 0, (1.7)where α = 2λ/(φ0β) and M = W0/C is gas injection Mach number. Theunstable mode happens where k → ∞ and where s is independent of k. Inthis cases ≈ φ0α2M2η(1− φ0M2)k2 (1.8)and corresponds to a growing perturbation with a weak oscillation. Thisinstability is associated with Bernoulli driving effect as growth rate dependson the square of gas injection speed M2, which in turn governs dynamicpressure variations in gas that scale as ρw2. Figure 1.6 shows oscillations indisplacement of the top of the magma column for various values of M, λ andη. The oscillations are the superposition of two modes: i. initial dominanthigh frequency oscillation and ii. most unstable long period mode. For mod-erate Mach number M = 0.1, the high-frequency oscillations are dominantfor a period of time and eventually exponential growth of the low-frequencyoscillations becomes significant. But as M is increased, the low frequencyoscillations grow faster before many of the high-frequency oscillations can13Figure 1.6: From [6] with permission. Oscillations in displacement uof the magma column at the column top, as predicted by thelinear stability analysis. The oscillations are the superpositionof only two modes, the fundamental mode and the least stablemode. Figures shown are for various values of M, λ and η asindicated. All frames have the same values of φ0 = 0.7, β = 100and γ = 0.occur.Non-linear wagging stability analysis is investigated in [6] for various pa-rameters, e.g. different gas flux injections and Mach numbers, by the sameapproach as linear analysis. The study [6] leads to qualitatively similarresults as the linear stability analysis. They show depending on the pa-rameters, the oscillation damping can be slow or moderately fast. As anexample, Fig. 1.7 shows the top displacement of the column for different14Figure 1.7: From [6] with permission. Oscillations in displacement uof the magma column at the column top, as predicted by thenon-linear stability analysis. The oscillations are the superpo-sition of only two modes, the fundamental mode and the leaststable mode. Figures shown are for various values of M, λ andη as indicated. Both frames have the same values of β = 100and γ = 2× 10−5.values of η, η = 5, 50. Higher η results in faster decay of high-frequencyoscillations.Observations of correlated tremor and gas flux variations show thatshort-period pulses in gas flux are nearly synchronous with tremors as audi-ble chugging [5, 22] and longer-period signals of average tremor amplitudeand gas flux are correlated with a nearly a minute lag [30] (see Figs. 1.8 and1.9). Bercovici et al. [6] show that short-period pulses are due to densityoscillations related to oscillations in column displacement since they providethe restoring force for the column wagging. Also, longer-period signals cor-respond to the Bernoulli effect of gas flux through the annulus at the end ofthe growing instability, which provides the driving mechanism to keep thetremor excited for long periods of time, see Fig. 1.6.Figure 1.8 also shows the growing sequence of beating envelopes oftremor. Bercovici et al. [6] suggest that the addition of harmonics of oscil-lation causes these envelopes. “Non-linear interactions between oscillations15afFigure 1.8: Seismic and Acoustic time series for Karymsky and San-gay volcanoes [22] show growing sequence of beating envelopesof tremor as well as correlated tremor and gas flux variations(acoustic time series).of a given frequency lead to cascading of energy to higher harmonics, andgiven the presence of both odd and even non-linearities, all harmonics (botheven and odd) are excited” [6]. The envelop structures are a superpositionof these close frequencies/harmonics, and the intervals between these en-velopes correspond to growth time for the unstable low frequency/Bernoullimode, see Fig. 1.10.When the magma column is displaced from its resting position, result-ing gas pressure variations within the annulus are imparted to conduit wallsand transmitted to seismometers. Two dimensional (linear) wagging is, thus,predicted to produce tremor signals on opposite sides of vent that are about1/2 cycle out of phase. It generates compressive ground motion at the di-rection of displacement and tensile ground motion on the opposite side. Thesign of ground motion around the volcanic vent will be reversed at half aperiod of oscillation later. Bercovici et al. [6] suggested that ground motionon opposite sides of a volcano should be spatially and temporally correlated.16(𝑎𝑎)(𝑏𝑏)Figure 1.9: From [30] with permission. a. Time series for both 10-saveraged tremor amplitude, and SO2 vent gas emission rate (inkg/s) for Fuego volcano in 2009. The top and bottom abscissaare offset by 32 s to compensate for overall lag. b. Evolution ofthe time lag between tremor and gas flux, where larger diamondsindicate correlation between time series larger than 0.65.Cross-correlations of seismic signals measured at two stations across the Re-doubt vent prior to its 2009 March eruption and showed that tremor signalsare out of phase by approximately half the dominant oscillation period (seeFig. 1.11).1.3.4 3D magma wagging modelBuilding on [6, 21], which are restricted to two-dimensional lateral waggingin Cartesian geometry, Liao et al. [26] extend the model to three dimen-sions, and column displacement can occur in any direction and plane, seeFig. 1.12. They identify new “whirling” regimes in the limit of very longvertical wavelength, i.e., k → 0 and viscous damping is negligible [21]. In17(𝑎𝑎)(𝑏𝑏)Figure 1.10: From [6] with permission. Time series and power spec-trum for surface displacement for a case with λ = 0.5, M = 0.1,and η = 1. The power spectrum is from a Fourier analysis ofthe time series with the background growth removed; power isin terms of U2(ω) where U(ω) is the discrete Fourier transformof U(H), and ω is the discrete angular frequency.this condition, displacement is continuous along the z direction, and at anyarbitrary height, the horizontal section of the magma column orbits, circularor elliptical orbit, around the center of the conduit. The elliptical/circularorbit is an expression of conservation of specific energy E0, and the shapeof the orbit is determined by the ratio of the linear component of specificenergy to the rotational component of specific energy ER (see §2.5). Oncethese energy components are determined, the solution can be found.Before an eruption, the increasing rate of gas exclusion and driving forcefrom the bubbles contribute to the total specific energy E0 accumulation,18Time lag (s)CorrelationFigure 1.11: Unbiased cross-correlation function between signals fromseismometers located north (RDN) and south (RSO) of Re-doubt volcano prior to its eruption on 2009 March 22. Thereis a positive correlation between the signals with a time lag of0.55 s. The time lag is approximately half the period of theperiod of the dominant signal (1 s). Dash lines show the 95%confidence interval for cross-correlation calculations.which can modify linear and rotational components of energy and, in turn,the wagging orbits. Figure 2.11 shows different elliptical whirling orbits.When the ratio of linear kinetic energy to rotational kinetic energy is zeroER = 0, whirling happens with a circular orbit. Increasing ER changeswagging toward more elliptical orbits. At very high energy ratios, ER→∞,whirling orbit changes to a straight line and the two-dimensional magmawagging in a single plane is resolved.Similar to two-dimensional magma wagging (where ER→∞), the radia-tion pattern of a seismic signal emanating from a three-dimensional magmawagging is correlated with motions of the magma column (see Fig. 1.13).Indeed, Liao et al. [26] suggested a method to identify the behavior of themagma column using the seismic data received from a series of seismic sta-tions around a volcano. The approach is based on calculating the relativetime-lag between seismic station pairs.19Figure 1.12: From [26] with permission. Sketch of the 3D magmawagging model. The column is assumed to have a free surfaceat the top and is coupled viscously to the magma system below.φ is the polar angle of the magma columns displacement, uis the radius displacement, or equivalently the displacementmagnitude. Rc and Rm are the radii of the magma columnand volcanic conduit, respectively.1.3.5 Summary of magma wagging model predictionsOn the basis of previous studies on magma wagging model, which is pre-sented briefly in the previous section §1.3.1, in this section we list the pre-dictions of magma wagging model that are testable and can be measured indata from field studies:1. The fundamental frequency of oscillation of a magma column inside afoamy vesicular annulus corresponds to volcanic tremor frequency (seeFig. 1.4 and equation 1.2).2. Increasing the annulus thickness causes the wagging frequency to de-crease as ω0 ∝ (R2c −R2m)−1/2 (see Fig 1.4).20𝑎𝑎)𝑐𝑐)𝑏𝑏)𝑑𝑑)Figure 1.13: From [26] with permission. a. Contours of phase-shift of P-waves generated by a 2D waggingat a frequency of 3.3 Hz. Black lines and dashed lines are two groups of equal-phase contours whichare out-of-phase with each other. Pentagons indicate the location of five virtual seismic stations.The orange dash line with arrowhead indicates the direction of the wagging plane.b. Contours ofphase-shift of P-waves generated by a counter-clockwise wagging with circular whirling orbit at afrequency of 3.3 Hz. Black solid lines and blue dashed lines are two groups of equal-phase contourswhich are out-of-phase with each other. Pentagons indicate the location of five virtual seismicstations. The orange dash line with arrowhead indicates the counter-clockwise wagging direction.c. P-waves calculated at the five virtual stations in a, which are in-phase or out-of-phase with eachother. Stations that are not in the same phase have half a period of oscillation T0/2 = 0.15 s lagrelative to each other.d. P-waves calculated at the five virtual stations in b.213. Increasing bubble volume fraction in the annulus causes the waggingfrequency to decrease as ω0 ∝√1/φ0.4. Toward the onset of volcanism when fragmentation disturbs annulusand porosity φ0 in the annulus vary in space, narrow band tremorsignal changes to broadband tremor as ω0 ∝√1/φ0, where φ0 variesin space.5. Degassing and audible chugging and tremor signals should be corre-lated (Figs. 1.6 and 1.8).6. Envelopes in tremor signals correspond to the superposition of modesthat are excited due cascading of energy from wagging frequency tothe higher harmonics (Figs. 1.10 and 1.8).7. Ground motions and related pressure waves distributed around thevolcano as the results of magma column displacements inside the vol-cano should be spatially and temporarily correlated (Figs. 1.11 and1.13).8. En route to the eruption magma column whirling orbit change becausemagma column energy changes, i.e., change in degassing behaviorand gas flux in the annulus changes the energy components of magmacolumn (Fig. 2.11).1.4 Critical knowledge gapsJellinek & Bercovici [21] propose the magma wagging model and explore itsproperties in 2D while the gas flow is not allowed in the annulus. Later,Bercovici et al. [6] extend magma wagging studies by allowing gas to flowinto permeable annulus and show how Bernoullis effect of gas flux can drivemagma waging. These two studies restrict the motion of the magma columnto a two-dimensional in a plane. Building on the work, Liao et al. [26] pro-pose an extended magma wagging model in 3D with impermeable annulusthat relaxes the 2D constraint and explore the spatial and temporal featuresof the seismic radiation pattern associated with tremor.22These studies are able to explain the ubiquity, persistence, and temporalbehavior of pre- and syn-eruptive volcanic tremors observed at different vol-canoes. Magma wagging theory is, thus, a potentially fully predictive modelfor explosive volcanic eruptions. A key knowledge gap, however, is exper-imental verification that the main underlying physical processes identifiedtheoretically are an accurate and complete description of the phenomenon.Critically, whether the long period Bernoulli mode can continuously drivethe predicted magma wagging is unclear. In this theory, the delivery of en-ergy from long period Bernoulli modes to short period wagging is only hy-pothesized. This energy transfer to do the work of wagging against viscousdissipation must be shown before the model can be taken to be reliable. Fur-thermore, additional processional and possibly torsional modes may emergetogether with the lateral wagging modes that are identified in 3D cylindricalgeometries. Perhaps most important related to explosion and fragmentationis that porosity can be modified and lead to non-uniform gas flux aroundthe magma column. With the aid of experiments, we can investigate how aninherent asymmetric distribution of gas flux modulates the Bernoulli modeand affects the wagging model, in turn.1.5 Thesis objectives and outlinesThe primary goal of this thesis is to test the magma wagging model exper-imentally. Next, we extend the model to explore spatially non-uniform gasflux, which is not tractable analytically and very challenging in 3D numer-ical simulations. To this end, we design and build an analog that capturesthe essential dynamics expressed through equations 1.1-1.6.The outline of this work is as follows: in Chapter 2, we first introduce theexperimental setup. Next, we discuss the mathematical modeling used foranalyzing the experimental results. We finish the chapter with a discussionof the dynamic similarity between experiment and the model applied to realvolcanoes and discuss the scaling in the problem. In Chapter 3, we presentexperimental results. We identify three regimes for the experiments andanalyze the mechanical conditions favoring each class of wagging. We end23this chapter with a discussion of the experimental results. In Chapter 4, tobetter understand and extend aspects of the experimental results, we carryout 2D Cartesian numerical simulations of our experiments. In Chapter 5we discuss the implication of our work for real volcanic activities. Finally,we conclude in Chapter 6 with a summary of the main results from thisthesis, further discussion of limitations of this study, and future work.24Chapter 2MethodsIn this chapter, we present the experimental setup and strategy. We thenintroduce some mathematical modeling we will use in this work to char-acterize experimental data. Finally, following this quantitative analysis ofour experiment, we introduce dimensionless control parameters that we useto understand our results as well as their relevance to establish dynamicsimilarity with the magma wagging model.2.1 Experimental setupOur analog experiments are carried out in a 10.16 cm (4′′) diameter× 180 cmlong hollow acrylic cylinder, which simulates a volcanic conduit. We install aviscoelastic rubber column with physical properties, lengths, and diametersthat we vary at the center of a manifold that forms the tank base. Whereasthe top of the column is a free surface, the column is fixed at the base witha peg that prevents vertical, radial and azimuthal motions. The rubbercolumns are made of Smooth-On’s ”VytaFlexTM” urethane rubbers in 50Aand 60A shore hardness with physical properties given in table 2.1. In theseexperiments the balance between restoring gas spring and viscous dampingforce expressed through η equation 1.7 is captured through prescribed phys-ical properties in our rubber column. This approach enables us to explorethe essential physics of magma wagging at the laboratory scale. It is im-25Column C1 C2 C3 C4D (cm) 8.9 8.9 7.6 7.6L (cm) 108 80 112 90ρ (kg/m3) 1040 1040 1040 1040Shore hardness 60A 50A 50A 60Aµ? (Pa.s) 8.18e4 2.02e4 1.59e4 1.22e4E? (MPa) 1.13 0.50 0.23 0.10Table 2.1: Physical properties of VytaFlexTM columns used in theexperiments. ? comes from experiments, see §2.3.2possible to capture both the elasticity of a permeable annulus and columninertia at the laboratory scale. Similarly, through this setup we capture thecritical damping effects that enable us to explore whether a Bernoulli modecan sustain higher frequency wagging.At the start of each experiment, we insert air at the base of the cylinderthrough the manifold with 12 equally-spaced 0.95 cm (3/8′′) diameter holespositioned around an 8 cm circle. We explore the responses of columns withprescribed elastic and linear damping properties to variously forced annularairflows. We control the distribution of air and the flow rate from each holeusing Panel-Mount flow-meters to measure the inflow. These flow-meterscan quantify the flow in the range of 4.2− 42 SCFM (1.982− 19.82 lit/sec)with 2% accuracy.To characterize and understand quantitatively, the excitation, evolution,and steady-state oscillating behaviors of our columns over a broad range ofconditions, we use three CANON EOS 80D cameras mounted normal toand above the columns, see Fig. 2.1. Each of these cameras records columnmovement with a rate of 29.97 frames per second. We put the hollow cylinderinside a rectangular tank and fill the volume between them with water. Thisapproach removes the challenge imaging through a curved boundary as waterand Acrylic glass have similar refractive indexes. Besides, we attach one Tri-Axis accelerometer on each side of the tank to capture the pressure variationinduced by column movement during each experiment. We use the Arduino26computing platform to control and record the accelerometers. The samplingrate of the data that these devices provide is around 90 Hz, see Fig. 2.1.2.2 Column calibrationWe carry out initial experiments to characterize the elastic and dampingproperties of each column. We also test the apparatus to ensure steady-state oscillation can be obtained. In doing this, we identify three columnwagging regimes during each experiment: i. a transient startup with theonset of gas forcing, ii. a statistical steady-state; iii. and a transient decaywhen the forcing is removed.We define the initial transient regime as the period from when the flow isintroduced around the rubber column until a steady response of the columnis established. The steady-state regime is marked by statistically station-ary column behavior and is the period between the end of the transientregime until the cessation of forcing. Oscillations in the steady-state regimeconstrain the inertial and elastic forces and reflect a balance between powerdelivered by gas flux and column dissipation. The damping regime followingcessation of forcing defines the viscous response of the column (Fig. 2.3). Weuse oscillation and displacement behavior in the damping regime to charac-terize the elastic and viscous properties of the rubber columns as they are ex-pressed through (e-fold) damping time scale, Td, for motions to cease and thenatural frequency of oscillation, ωn (§2.3.1). For example, Fig. 2.2 a. showsan evolution of a trajectory in x − y plane of top of the column, which iscaptured by the camera installed above the experimental setup, normalizedby the gap width between the hollow cylinder and rubber column δ in x andy direction. The data show a short initial transient, a sustained period ofsteady-state behavior and a transient decay in amplitude where gas forcingis turned off.In Fig. 2.3 we highlight the initial transient and damping periods. Fig-ure 2.3 a. and d. are examples of the trajectory of top of the column in thetransient and damping regime, respectively. When input flow is introducedaround the column, it displaces from the resting position, and the displace-27Hollow cylinderSurrounding boxFrameAccelerometerRubber ColumnTop view cameraSide view camerasddddd 𝑥𝑥𝑦𝑦𝑜𝑜𝑂𝑂1𝑟𝑟𝑂𝑂12𝑂𝑂11𝑂𝑂10𝑂𝑂9𝑂𝑂8𝑂𝑂7𝑂𝑂6𝑂𝑂5𝑂𝑂4𝑂𝑂3𝑂𝑂2dFigure 2.1: Experimental setup and the arrangement of inlets inthe manifold. Hollow cylinder dimensions: 10.16 cm inter-nal diameter × 180 cm long. Surrounding box dimensions:20 cm×20 cm×125 cm. The top camera is fixed 20 cm above theapparatus, and the distance between side cameras and the tanksurface is 150 cm. Distance between the center of the manifoldand inlets r is 4 cm, and inlets diameter d is 0.95 cm.28(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)Steady-State DampingTransient𝛿𝛿(𝑑𝑑)Figure 2.2: a. Example of trajectory of top of the column. Displacements in x and y directions are scaledby δ. Colorbar shows the dimensionless time t/Td. b. The corresponding time series in x direction,and c. The corresponding time series in y direction. Transient-startup, steady-state, and dampingregimes are defined in the text.29ሺ𝑎ሻ ሺ𝑑ሻሺ𝑏ሻሺ𝑐ሻ ሺ𝑓ሻሺ𝑒ሻFigure 2.3: a. Example of trajectory of top of the column duringthe transient forcing, b. the corresponding time series in xdirection, and c. the corresponding time series in y direction. d.Example of trajectory of top of the column during the dampingregime, e. the corresponding time series in x direction, and f.the corresponding time series in y direction.ment amplitude increases until the steady-state regime emerges. During thedamping regime (Fig. 2.3), after airflow stops, oscillations decay in ampli-tude and column generally recovers a resting position. Therefore, it oscillateswhile the amplitude of oscillation decreases exponentially over time until itreaches the resting position. Corresponding time series to the transient anddamping trajectories are shown in Fig. 2.3 b-c and e-f, respectively.30Figure 2.4: Schematic of simplified version of experiments viewedfrom top.2.3 Calibrating the column response propertieswith mathematical models2.3.1 2D mass-spring-damper modelFrom Fig. 2.2 and Fig. 2.3, the column has a characteristic damping time,Td, and an oscillation frequency, ωn, that are consistent during the initial,steady-state, and damping regimes. Thus, here we develop mathematicalmodels for the response of the system to the imposed gas forcing to constrainthe elastic properties of the column quantitatively. This exercise enable usto verify whether our experiments adequately captures the effect of η inequation 1.6.When we look at the column from the top, we can simulate the systemto a two degree of freedom mass-spring-damper system, see Fig. 2.4. In thissystem c is damping coefficient, K is spring constant, and M is the mass inthe system. Newton’s second law of column motion can be written,d2udt2+cMdudt+KMu = f, (2.1)31where u = u(x, y) is the displacement and f is the effective force per mass.equation 2.1 is usually written into the form of:d2udt2+ 2ζωndudt+ ω2nu = f, (2.2)where the damping ratio ζ and natural frequency ωn of the system are:ζ =c2√MKand ωn =√KM.For a free oscillation akin to the damping regime (f = 0), the solution ofequation 2.2 is:u = Ae−ζωntcos(ωdt+ φ), (2.3)where the amplitude A and phase φ determine the behavior needed to matchinitial conditions, and the fundamental damped frequency ωd is:ωd = ωn(1− ζ2)1/2. (2.4)We non-dimensionalize equation 2.2 by selecting the gap width δ as thelength scale and the column relaxation time Td = 1/ζωn as the dampingtime scale. In this case we define u(x, y) = u′(x, y)δ, t = t′Td and theequation of motion changes to (after dropping the primes):d2udt2+ 2dudt+ γ2u = f, (2.5)where the balance between elastic and restoring damping forces (akin to ηin equation 1.6) is expressed through the ratio of time scales for dampingand natural oscillation:γ = Tdωn. (2.6)We solve equation 2.5 during transient decaying regime where f = 0 byassuming u(t) = ert and it results in:u(t) = C1e(−1+i(√γ2−1+φ))t + C2e(−1−i(√γ2−1+φ)t, (2.7)32where C1, C2, and φ are determined by imposing the initial conditions. Itis to be observed from equation 2.7 that decaying constant is unity and thevalue of γ determines the behavior of system in the free oscillation:• Overdamped γ < 1: The system exponentially decays to steady-statewithout oscillating.• Critically damped γ = 1: The system returns to steady-state asquickly as possible without oscillating.• Underdamped γ > 1: The system oscillates with a slightly differentfrequency than the natural frequency with the amplitude graduallydecreasing to zero.As we discussed in §2.2, we characterize the damping properties of rub-ber columns using the transient decaying regime. For more detail, assumingthat the column damping is linear (equation 2.5 holds), the amplitude ofoscillation will decay as exp(−t/Td). From experimental fits to the displace-ment time series in the damping regime, we use an exponential fork envelopof the oscillation to find the decay rate of the amplitude in the free oscillation1/Td and the characteristic damping time.Figure 2.5 a-b are examples of fitted curves to damping regime in x and ydirection for one of the experiments. Calculated damping coefficients (ζωn)is 0.20 in x and y direction. Figure 2.6 shows the distribution of the dampingcoefficient calculated from all the fitted curves to damped regimes duringexperiments. Solid lines are fitted Gaussian curve to the distributions. Themean values of damping coefficients ζωn and its standard deviation σ, fordifferent columns can be read from this figure and is reported in Table 2.2.From spectral analysis of time series of the damping regime we also find,in turn, the fundamental damped frequency (ωd = 2pifd). Figures 2.7 a. andb. show the frequency spectrum of damping time series in Figs. 2.5 a. andb., respectively. Both time series demonstrate peak frequency at 0.29 Hz.Figure 2.8 shows the distribution of fundamental damped frequency cal-culated from all the damped time series during experiments. Solid lines arefitted Gaussian curve to the distributions. The values of damped frequencies33(𝑎) (𝑏)𝑒−𝜁𝜔𝑛𝑡+𝑏 𝑒−𝜁𝜔𝑛𝑡+𝑏−𝑒−𝜁𝜔𝑛𝑡+𝑏−𝑒−𝜁𝜔𝑛𝑡+𝑏Figure 2.5: a. The trajectory of top of the column in x direction, b.The trajectory of top the column in y direction. Dashed linespresent envelopes fitted to trajectories to calculate damping co-efficient.Column ζωn σC1 0.42 0.2C2 0.38 0.15C3 0.15 0.03C4 0.2 0.19Table 2.2: Values of damping coefficients and its standard deviationfor different columns.(ω¯d) and its standard deviation (σ) for different columns can be read fromthis figure and is reported in Table 2.3.With the linear damping we can calculate the natural frequency of freeColumn ω¯d σC1 3.38 0.61C2 4.33 0.17C3 2.10 0.17C4 1.62 0.20Table 2.3: Values of fundamental damped frequency and its standarddeviation for different columns.34(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)(𝑑𝑑)Figure 2.6: Histograms of damping coefficient (ζωn in 2.3) calculatedfrom damped time series. Solid lines represents the fitted Gaus-sian curves. a., b., c., and d. correspond to C1, C2, C3, andC4 respectively.oscillation ωn by substituting Td = 1/ζωn in equation 2.4 and rearrangingit in the form of:ω2n = ω2d +1T 2d. (2.8)We have used two-dimensional mass-spring-damper model to character-ize relaxation time, Td, and oscillation frequency, ωn, of columns used inexperiments. The values of Td, ωn, and γ for different rubber columns usedin our experiments are shown in Table 2.4. Note that the transient-startupregimes usually takes between 1 to 3 relaxation times, e.g., see Figs. 2.2 and2.3 but we find statistically stationary steady-state is related to experiments35(𝑎𝑎) (𝑏𝑏)Figure 2.7: a. power spectrum of Fig. 2.5 a., b. power spectrum ofFig. 2.5 b.(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)(𝑑𝑑)Figure 2.8: Histograms of damped angular frequency (ωd) calculatedfrom damped time series. Solid lines represent the fitted Gaus-sian curves. a., b., c., and d. correspond to C1, C2, C3, andC4 respectively.36Column T¯d ω¯n γC1 2.38 3.40 8.16C2 2.63 4.35 11.53C3 6.66 2.10 14.08C4 5.00 1.63 8.15Table 2.4: Values of relaxation time, natural frequency of the system,and column’s characteristics number.𝑢𝑢𝑦𝑦𝑧𝑧𝑥𝑥𝑑𝑑𝑧𝑧𝑁𝑁𝑁𝑁𝑑𝑑𝑁𝑁𝑅𝑅𝑁𝑁𝑁𝑁𝑟𝑟𝜃𝜃Figure 2.9: Schematic of simplified version of experiments from sideview.after 5 relaxation times Td, in particular.2.3.2 Cantilever beam modelWhereas equation 2.5 provides insight into relaxation time and natural fre-quency of the column, we need to study column’s bending to understandthe elastic properties of column and to constrain the effective viscosity andYoung’s modulus of the column.When we look at the column from side, we can simulate it to a viscoelas-tic cantilever beam. We solve cantilever beam problem to characterize the37viscous properties of the column and its elastic response to the forcing. Togovern the equation of motion for the column, we consider a viscoelasticbeam that has a cross-sectional area A, moment of inertia I, Young’s modu-lus of elasticity E, material density ρ, and efficient viscosity µ. Now consideran element of the length dx of the beam that is bent on one side, e.g. Fig. 2.9.Here, the column is compressed at the right side of neutral axis NA and isdecompressed at the left side of NA. Strain  at distance r from NA canbe written, =(R− r)θ −RθRθ=rR, (2.9)where R is the radius of bending curvature. We use KelvinVoigt modelconsists of a Newtonian damper and elastic spring connected in parallel tomodel viscoelasticity of the material. Stress σ related to strain  in thismodel can be written,σ = E+ µ˙, (2.10)where ˙ is the strain rate and from equation 2.9 it can be written,˙ = − rR2R˙. (2.11)From the definition of stress σ we have:df = σ dA, (2.12)where df is the force applied to the area dA, and the partial bending momentdM from this force is:dM = r df, (2.13)now we can find bending moment M by integrating dM across the area A:M =∫σrdA =∫(ER− µR2R˙)r2dA = (ER− µR2R˙)I. (2.14)Assuming small deflections along the column, bending curvature R can bewritten,1R=∂2u(z, t)∂z2, (2.15)38𝑑𝑑𝑧𝑧𝑉𝑉(𝑧𝑧, 𝑡𝑡)𝑀𝑀(𝑧𝑧, 𝑡𝑡)𝑉𝑉 𝑧𝑧, 𝑡𝑡 + 𝜕𝜕𝑉𝑉 𝑧𝑧, 𝑡𝑡𝜕𝜕𝑧𝑧𝑑𝑑𝑧𝑧𝑤𝑤(𝑧𝑧, 𝑡𝑡)𝑂𝑂𝑀𝑀 𝑧𝑧, 𝑡𝑡 + 𝜕𝜕𝑀𝑀 𝑧𝑧, 𝑡𝑡𝜕𝜕𝑧𝑧𝑑𝑑𝑧𝑧Figure 2.10: Free body diagram for an element of column. w is ex-ternally distributed load, V is shear force, and M is bendingmoment.Where u(z, t) is column deflection at distance z and time t, hence,− R˙R2=∂3u(z, t)∂t∂2z, (2.16)where t represents time. Then substituting equation 2.15 and equation 2.16in equation 2.14 gives the relation between bending moment and deflectionin viscoelastic beams,M(z, t) = EI∂2u(z, t)∂z2+ µI∂3u(z, t)∂t∂z2. (2.17)A free body diagram for an element of column is shown in Fig. 2.10. Ifwe derive Newton’s second law for this element we have:ρAdz∂2u(z, t)∂t2= w(z, t)dz +∂V (z, t)∂zdz, (2.18)then for the element, neglecting rotary inertia, taking moments about O:∂M(z, t)∂z= −V (z, t), (2.19)∂V (z, t)∂z= −∂2M(z, t)∂z2, (2.20)39finally substituting equation 2.17 and equation 2.20 into equation 2.18 leadsto the governing equation of motion for a viscoelastic cantilever beam undera distributed load w(z):∂4u(z, t)∂z4+µE∂∂t∂4u(z, t)∂z4+ρAEI∂2u(z, t)∂t2=w(z, t)EI. (2.21)We non-dimensionalize equation 2.21 by choosing the gap width δ as thelength scale for y direction, column height L as the length scale for z direc-tion, and Td as the time scale. We define u = u′δ, z = z′L, t = t′Td and theequation of motion changes to (after dropping the primes):∂4u(z, t)∂z4+µETd∂∂t∂4u(z, t)∂z4+ λ4∂2u(z, t)∂t2= w(z, t), (2.22)where,λ4 =ρAL4EIT 2d(2.23)and λ is the characteristic of the column.We are now going to solve equation 2.22 during transient decaying regimewhen w(z, t) = 0 by multiplying it by eikz and integrating with respect toz. Taking the Fourier transform with respect to z, for each fixed t, of u(z, t)by:uˆ(k, t) =∫ +∞−∞u(z, t)e−ikzdz, (2.24)and applying this solution to equation 2.22,λ4∂∂t2uˆ(k, t) + k4µETd∂∂tuˆ(k, t) + k4uˆ(k, t) = 0, (2.25)we can solve equation 2.25 by trying zˆ(k, t) = ert, which r can be a complexnumber, and it results in,uˆ(k, t) = e−αte±iωdt, (2.26)40whereα =k42λ4µETdand ω2d = 2αETdµ− α2,so the general solution to equation 2.25 is:uˆ(k, t) = Uˆ1(k)e−αte+iωdt + Uˆ2(k)e−αte−iωdt, (2.27)where Uˆ1(k) and Uˆ2(k) are arbitrary constants. Taking the inverse Fouriertransform of equation 2.27 we obtain a linear ODE for u(z):∂4u∂z4− 2αETdµλ4u = 0, (2.28)and the general solution to the beam equation equation 2.28 is:u = C1 cos λ′z + C2 sin λ′z + C3 cosh λ′z + C4 sinh λ′z, (2.29)whereλ′4 =2αETdλ4µ, (2.30)and constants C1, C2, C3, and C4 are determined from boundary conditions.Applying the boundary condition for the cantilever beam:u = 0,dudz= 0 at z = 0,d2udz2= 0,d3udz3= 0 at z = 1,non-trivial solutions are found to exist only if cosh(λ′) cos(λ′) + 1 = 0 or:λ′ = 1.875, 4.694, 7.885, ... .which different λ′ correspond to different modes of oscillation.As we mentioned at the beginning of this section the ultimate goal hereis to characterize the viscoelastic properties of the column, and we knowthat time dependent solutions for equation 2.5 and equation 2.22 for thefirst mode of oscillation must be the same. By comparing equation 2.7 and41equation 2.27 we have:α = 1 and ω2d = γ2 − 1, (2.31)hence effective viscosity µ and Youngs’ modulus can be found by comparingequation 2.23, equation 2.30, and equation 2.31,µ = 2ρAL4ITd11.8754and E =µγ22Td, (2.32)where ρ is column density. Calculated values of E and µ for different columnsare presented in Table 2.1. Note that if the column was pure elastic µ = 0,decay constant vanishes α = 0, the oscillation frequency becomes the naturalfrequency ω2n =λ2EIρAL4and equation 2.27 changes to:∂4u(z, t)∂x4− λ4u(z, t) = 0, (2.33)which can be solved as equation 2.28.2.4 Scaling considerationsFor each experiment, different inflows are introduced from each of the 12inlets to the cylinder. Columns respond differently with respect to the totalinflow and the arrangement of inflows. To characterize the symmetry andintensity properties of these different forcing scenarios and the mechanicalresponse of the column, we define 2 parameters: a forcing eccentricity e andDeborah number De. We usee =r′r, (2.34)as a metric for the asymmetry of the flow around the column. Here, r isthe distance between the manifold center and the inlet center, r = |Oi−O|,Fig. 2.1. The position r′ is the distance between weighted average positionof all inlets and the center of the manifold:r′ =√x¯2c + y¯2c ,42wherex¯c =Σ12i=1QixiΣ12i=1Qiand y¯c =Σ12i=1QiyiΣ12i=1Qi.Here (xi, yi) and Qi are the ith inlet position and its flow rate. On thebasis of this definition, e changes between 0 and 1, 0 presents the perfectsymmetric forcing and 1 is the most asymmetric forcing, e.g. there is onlyone active inlet.To characterize the amplitude of forcing for a specified eccentricity e andcolumn damping we use a Deborah number:De =TdTe, (2.35)where Td is the viscoelastic relaxation time of the rubber column and Te isthe time scale for the build up of elastic stresses in response to the forcing.Balancing the power delivered to an oscillation as a result of an imposedflow at mean velocity v¯ with the build up of elastic potential energy givesus:∆KE∆t∼ PETe, (2.36)where KE = ρ0v¯2/2 is the flow kinetic energy density, ∆t is a timescale fordelivery of this energy to do the work of orbiting an oscillation, and PE isthe column elastic potential energy density related to the build up of strain(δ/L). The average velocity v¯ of the flow around the column is given by:v¯ =Σ12i=1QiA, (2.37)where A is the area between rubber column and hollow cylinder:A =pi4(D2 − (D − δ)2), (2.38)and D is the internal diameter of hollow cylinder. If we take the advectivetime as ∆t = v/L then:∆KE∆t=12ρ0v3L. (2.39)43To find PE we assume that column is under a constant distributed loadsuch as maximum deflection on top of the column is δ. The correspondingload w will be:w =8EIδL4, (2.40)bending moment M along the column due to distributed load w is:M =w(L− z)22, (2.41)where z is the distance from the fixed side. Then bending energy U storedin the column can be calculated from:U =∫ L0M22EIdz =8EIδ25L3. (2.42)Now we can find the column elastic potential energy density PE by:PE =UV, (2.43)where V is the volume of the column. Combining equation 2.39, equa-tion 2.42, and equation 2.43, the elastic time scale is,Te =645piEρ0Iδ2L3D2v¯3. (2.44)equation 2.44 shows that elastic response time of the column Te is propor-tional to its Young’s modulus of elasticity E and represents the importanceof spring behavior of the column. On the other hand, relaxation time Td isproportional to effective viscosity of the column µ and shows the importanceof viscous damping. As a result, Deborah number De = Td/Te indicates theimportance of viscous forces in the column to elastic behavior of the col-umn which is similar to the key dimensionless number η = νm/(CU) in thewagging model, equation 1.6, which is a metric for importance of viscousdamping in the magma column to the gas spring force in the annulus inthe original theoretical magma wagging model. In previous studies [6, 21] ηvaries in range of 5 − 50. In this study, Deborah number changes between4410 − 3 × 104. There is a notable difference between ranges of De and ηbecause of their different definitions. De is a function of column propertiesand gas velocity (De ∝ v¯3), however, η is defined by column properties andconstant sound velocity in the gas.2.5 Energetic contributions to column toptrajectoriesIn Chapter 3, we will show that the trajectory of top of the column variessignificantly for different airflow forcing. Following Liao et al. [26], wecharacterize steady-state trajectories in terms of contributions from linearand rotational kinetic energy as well as elastic potential energy.. Here, weuse energetic components to define an energy ratio ER, which we applylater as a means to distinguish between different column behaviors. Atsteady-state where the power delivered through the injected air is balancedby dissipation as a result of linear damping, equation 2.22 reduces to,d2udt2= −ω2nu, (2.45)In cylindrical-polar coordinate, equation 2.45 can be recast as two coupledevolution equations, for the radial and azimuthal displacement, respectively[3]d2udt2− udφdtdφdt+ ω2nu = 0, (2.46)ud2φdt2+ 2dudtdφdt= 0, (2.47)We define specific energy E0,E0 = KELin +KERot + PE, (2.48)whereKELin =12(du(t)dt)245Figure 2.11: Different elliptical whirling orbits corresponding to dif-ferent ER = KELin/KERot values. Circular orbit correspondsto ER = 0. When ER→∞ orbit changes to a straight line.is linear kinetic energy,KERot =12u(t)2(dφ(t)dt)2is rotational kinetic energy, andPE =12ω2nu(t)2is potential energy. By definition at steady-state the total energy is con-served and the ratio of linear kinetic energy KELin to rotational kineticenergy KERot determines the shape of column trajectory.ER =KELinKERot, (2.49)where KELin and KERot are time-averaged for each experiment. Fig-ure 2.11 shows experiments of different elliptical whirling orbits. Pure ro-tational motion/circular orbit occurs when ER = 0. By increasing ER,orbits become more elliptical. And, two-dimensional wagging oscillation isrecovered in the limit of ER→∞.46Chapter 3Experimental Results andDiscussionA given column displacement depends on the intensity and asymmetry ofthe air forcing from below. Here we study the displacement properties ofeach column to imposed forcing eccentricity e and intensity De for roughly1000 individual experiments by using methods introduced in Chapter 2. Inthis chapter, we first characterize trajectories of the top of columns qualita-tively in the steady-state regime into three cases of response of the column,depending on e−De conditions; “rotational”, “chaotic”, and “mixed-mode”.We then use the energy ratio ER to define the rationale for each response.Next, we examine the time series of the top trajectory spectrally and in-vestigate the correlation and coherence between time series of accelerationcaptured by accelerometers for each response. Finally, we show all 1000 ex-periments on one regime diagram to investigate the effects of intensity andasymmetry of the airflow to the column responses.3.1 Rotational, mixed-mode, and chaotic columnresponsesHere, we introduce “rotational”, “chaotic”, and “mixed-mode” responses forcolumn trajectories and displacement in the steady-state regime.47Figure 3.1 is an example of a rotational motion that emerges for a sym-metric airflow condition around the column with eccentricity e = 0 and aforcing intensity of De = 1.4 × 104. The top of the column in rotationalresponse follows a circular orbit with an approximately constant radius, seeFigs. 3.1 a, c, d., and Fig. 3.2. In addition to this motion of the column top,continuous vertical deflections over the column length occur during eachorbit (Fig. 3.1 b).By contrast, Fig. 3.3 is an example of mixed-mode motion. In thisexample, the airflow distribution around the column is slightly asymmetriccomparing to Fig. 3.1 and Fig. 3.2 with forcing eccentricity e = 0.08 andforcing intensityDe = 6.8×103. Under these e−De conditions the trajectoryof the column top orbits changes continuously over time (Fig. 3.3) with shiftsin direction. The column’s top follows predominantly circular trajectorieswith radial explosions over time scales < Td as well as time scales Td(Fig. 3.4). Averaged over the full length of the experiments, column motionsare, time varying combinations of rotational and linear displacements. Thecolumn intermittently ends up in rotational trajectories, but these responsesare transient and variable in their longevity. In addition, vertical deflectionsof the column over its length during each orbit can be continuous or periodic.Finally, the amplitude of oscillations in x/δ − y/δ space are not constantand characterized by emergent beating envelopes, Figs. 3.3 c. and d.Figure 3.5 is an example of a chaotic displacement response. In thisexample the airflow around the column is highly asymmetric comparing toFig. 3.1 and Fig. 3.3 with e = 1 and De = 75. Trajectories shift continuouslyover time scales < Td, ∼ Td, and Td with mainly linear displacementsthat are stochastic in time (Figs. 3.5 & Fig. 3.6). Vertical deflections of thecolumn are chaotic and quasi-periodic in space and similarly stochastic intime. Similar to Fig. 3.3 the amplitude of oscillation in x/δ− y/δ space arenot constant, generally smaller than in rotational and mixed-mode responseswith beating envelopes that are complex in their shape, where they exist atall.48(𝑎𝑎)(𝑐𝑐)(𝑑𝑑)(𝑏𝑏)01𝑧𝑧/𝐿𝐿𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕Figure 3.1: Rotational response; a. The trajectory of top of the column. Displacements in x and ydirections are scaled by δ. The colorbar is the dimensionless time t/Td. b. Snapshots of columndeflection during an orbit. c. The corresponding time series in x direction, and d. The correspondingtime series in y direction. Blue lines are upper and lower peak envelops of displacement time series.49Figure 3.2: Trajectory of top of the column in a rotational response. Here flow is evenly distributed aroundthe column e = 0 and Deborah number is De = 1.4× 104. Each subplot shows the trajectory duringone dimensionless time Td, and the colorbar shows the corresponding time. The column followsapproximately constant orbits 2 to 3 times in each time window. Trajectories are similar to eachother in each subplot, and the ones in other subplots in terms of radius and aspect ratio. Thedirection of rotation is constant, i.e. it moves clockwise.50(𝑎𝑎)(𝑐𝑐)(𝑑𝑑)(𝑏𝑏)01𝑧𝑧/𝐿𝐿𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕Figure 3.3: Mixed-mode response; a. The trajectory of top of the column. Displacements in x and ydirections are scaled by δ. The colorbar is the dimensionless time t/Td. b. Snapshots of columndeflection during an orbit. c. The corresponding time series in x direction, and d. The correspondingtime series in y direction. Blue lines are upper and lower peak envelops of displacement time series.51Figure 3.4: Trajectory of top of the column in a mixed-mode response. Flow distribution around the columnis approximately even e = 0.08 and Deborah number is De = 6.8 × 103. Each subplot shows thetrajectory during one dimensionless time Td and colorbar shows the corresponding time. The columnmainly follows elliptical and circular orbits, e.g. t = 25, 26, 33 are elliptical orbits while t = 34 and36 are circular. Aspect ratios and radius of ellipses and circles are different from time to time, e.g. theevolution of trajectory from t = 30 to t = 33 as the orbits are enlarging. Sometimes we observe linearmovements, e.g. t = 28. Although trajectories change over time but the direction of motion is alwaysthe same, i.e. for this example column moves in clockwise direction.52(𝑎𝑎)(𝑐𝑐)(𝑑𝑑)(𝑏𝑏)01𝑧𝑧/𝐿𝐿𝒕𝒕𝒕𝒕𝒕𝒕𝒕𝒕Figure 3.5: Chaotic response; a. The trajectory of top of the column. Displacements in x and y directionsare scaled by δ. The colorbar is the dimensionless time t/Td. b. Snapshots of column deflectionduring a cycle. c. The corresponding time series in x direction, and d. The corresponding time seriesin y direction. Blue lines are upper and lower peak envelops of displacement time series.53Figure 3.6: Trajectory of top of the column in a chaotic response. Flow distribution around the column ishighly uneven e = 1 and Deborah number is De = 75. Each subplot shows the trajectory during onedimensionless time Td and colorbar shows the corresponding time. Linear trajectories are the mostobserved pattern in this figure while the direction of linear trajectories are different, e.g. trajectoryat t = 41 is parallel to x axis while at t = 46 is parallel to y axis. At t = 40 and t = 44 we see thatdirection of linear motion changes. There are two elliptical trajectories at t = 37 and t = 39 whilethe orientation and size of ellipses are different. Other trajectories are a combination of linear andcircular movements in small scales as we see at t = 36 or t = 43.543.2 Energy distributionIn this section, we characterize steady-state column motions in terms ofthe time-varying partitioning of rotational, kinetic, and potential energy(definition in §2.5) in rotational, chaotic, and mixed-mode responses. Ourgoals are to: i. characterize each response in greater detail, ii. identify aquantitative metric for end-member rotational and chaotic motions as wellas the distinctive properties of the mixed-mode regime.Figure 3.7 shows time series of these three components of the specific en-ergy E0 for different motions. During a rotational response (Fig. 3.7 a., thetotal specific energy and its different components are statistically stationary,fluctuate around an approximately constant value, consistent with predomi-nantly rotational trajectories (Fig. 3.1), KERot has the highest contributionto E0 while KELin is the smallest component of the specific energy.During a mixed-mode response, by contrast, the total specific energy, aswell as each energetic component, varies significantly over Td as well as overmuch longer time scales. In general, KELin is the smallest component of E0over time scales ∼ Td and Td while energy is partitioned approximatelyequally between PE and KERot. A feature of the mixed-mode response isthat the column intermittently enters protracted periods marked by rota-tional motions that are also statistically stationary. For example, over thetime interval 39 6 t 6 44 motions are predominantly rotational and KERotis maximized.Finally, all three components of energy are comparable in magnitudeduring a chaotic response (Fig. 3.7 c). A marked difference to the rotationaland mixed-mode responses is that the contribution of KELin can exceedKERot.We use histograms of the distribution of specific energy componentsfor all the experiments in rotational, mixed-mode, and chaotic responserespectively, to distinguish the time-averaged properties of each response,Figs. 3.8 a, b, c. During rotational responses, KERot is the dominant com-ponent of energy with an average of 64% of total specific energy in a rangeof 55% < KERot < 76%. The stored PE is 32% on average with a range of55(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)Transient rotational responsesFigure 3.7: a. Energy time series in a rotational response. Flow conditions are e = 0 and De = 1.4× 104.,b. Energy time series in a mixed-mode response. Flow conditions are e = 0.08 and De = 6.8×103., c.Energy time series in a chaotic response. Flow conditions are e = 1 and De = 75. Time is scaled withTd and rotational kinetic energy KERot, linear kinetic energy KELin, potential energy PE, and totalspecific energy E0 are normalized by the maximum of specific energy E0max during each experiment.56Response KELin σKELin PE σPE KERot σKERotRotational 0.04 0.02 0.32 0.06 0.64 0.06Mixed-mode 0.10 0.01 0.29 0.05 0.61 0.05Chaotic 0.32 0.07 0.18 0.09 0.50 0.07Table 3.1: The average and standard deviation of time-averaged en-ergy components for all experiments in each response.18% < PE < 45%. The contribution from linear kinetic energy KELin 4%on average with a range of 1% < KELin < 8% (Fig.3.8 a).In a mixed-mode response, similar to the rotational mode, energeticcomponents are distinct in mixed-mode response and KERot is the dominantcomponent of energy with an average of 61% in a range of 51% < KERot <68%. The stored PE is 29% on average with a range of 22% < PE < 40%.The linear kinetic energy KELin contributes 10% on average and varies in7% < KELin < 12% (Fig.3.8 b).In contrast to the rotational and mixed-mode responses, energetic com-ponents of a chaotic modes are comparable to each other. The rotationalkinetic energy KERot with an average of 50% varies in 21% < KERot <71%. The potential energy PE with the average of 18% changes between1% < PE < 68%. The kinetic linear energy KELin contributes 33% onaverage with a range of 10% < KELine < 47% (Fig.3.8 c). We presentthe average value and standard deviation for time-averaged components ofspecific energy during all experiments in each mode in Table 3.1.To summarize, we observe in §3.1 that circular orbits are the main pat-tern in rotational responses and linear displacement is the dominant motionin chaotic motions. In addition, consistent with these kinematic observationswe show that the contribute of linear kinetic energy to total motion increasesform rotational to chaotic responses while rotational components of energyare still the dominant energy components. To quantify this behavior wenow use ER, defined in §2.5, as a metric to distinguish between different re-sponses. Figure 3.8 d. shows energy ratio ER = KELin/KERot distributionfor three responses and dashed lines represents the fitted Gaussian curves57(𝑎𝑎) (𝑏𝑏)(𝑐𝑐) (𝑑𝑑)Figure 3.8: Distribution of time-averaged specific energy componentsfor all the experiments in a. rotational responses, b. mixed-mode responses, and c. chaotic response. d. Distribution ofenergy ratio ER for different motions. Dashed lines are fittedGaussian curves to ER distributions, detail in Table 3.2.to the ER distribution in each one. The energy ratio ER for rotationalmotions varies in 0.01 < ER < 0.12 with the average of 0.06. Mixed-moderesponse is observed with ER in the range of 0.12 < ER < 0.22 and itsmean value is 0.16. Chaotic modes occur over a broad range of energy ratioER > 0.20. The average of ER for chaotic motions in our experiments is0.71. The mean value and standard deviation corresponding to the fittedcurves is represented in Table 3.2.58Response ER σRotational 0.06 0.03Mixed-mode 0.16 0.02Chaotic 0.71 0.19Table 3.2: The average value of energy ratio ER and its standarddeviation for different responses.3.3 Spectral analysisThe time series of trajectories in §3.1 and energetic analysis in §3.2 showthat rotational as well as linear oscillations occur in response to the forcedairflow. Accordingly, in this section, we use spectral analysis of both dis-placement and accelerometer time series to characterize the periodic partsof the motions in all three responses.Figure 3.9 shows the normalized power spectrum of time series of tra-jectory and normalized power spectrum of time series of accelerometer datafor three examples in different responses in black and red color, respectively.We normalize the frequency axis with the damped natural frequency of thecolumn, fd, to make clear the damping natural frequency (see Fig. 2.7) andthe higher frequency oscillation related to the air forcing. The frequencyresponse characteristics of accelerometers preclude our ability to resolve fre-quencies as low as the damped frequency of columns fd. Consequently, tocapture the full column response, we combine both time series. We recordvideos of top displacement with a rate of 29.97 fps. Therefore, its Nyquistfrequency is 15 Hz, while accelerometers record 90 samples in a second,and its Nyquist frequency is 45 Hz. And they can record higher frequencyvibrations than what a video can capture.Figure 3.9 shows that periodic motions exist in all three examples. Thereare three peak frequencies that are roughly at 1× fd, as well as 10× fd, and25 − 30 × fd. Note, fd is the natural frequency of the system (Fig. 2.7).For all the responses, the second and third peaks have lower power than thefirst peak, but the ratio of second and third peak power to the first peak for59chaotic response is greater than mixed-mode and rotational response.Figure. 3.10 shows the normalized power spectrum of all the experimentsin each response, a. rotational, b. mixed-mode, and c. chaotic. In rotationalresponses, the first peak corresponding to fd occurs at 1−2 ×fd. On average,the second higher frequency peak related to the air forcing is approximately10% of the first peak power (Fig. 3.10 a). Similar to rotational mode, in amixed-mode response, the damped natural frequency is at 1− 2 × fd and,on average, the second peak power is approximately 15% of the first peakpower (Fig. 3.10 b). The damping natural frequency in chaotic responsesis broader compared to rotational and mixed-mode regimes, and the secondpeak power is,on average, is 30% of the first peak power which has a highervalue than the other two responses.3.4 Regime diagramWe introduced rotational, mixed-mode, and chaotic response to different air-flow conditions in §3.1. Then we explored differences between these motionsin terms of the energy content of top trajectories in §3.2. Here we investigatethe column response in terms of the trajectory of top of the column and itsenergy content to the airflow eccentricity e and intensity De.Figure 3.11 summarize the response characteristics of all the experimentscarried out with four columns in our study in terms of trajectory type andenergy ratio ER to the airflow eccentricity e and intensity De. For sym-metric forcing e = 0, as De is increased, initially chaotic mode is dominant,transition to mixed-mode response emerge at De > 6×103, and the columnresponse evolve to be predominantly rotational for De > 104. If forcingeccentricity is increased to 0.5, e < 0.5, the transition from chaotic responseto mixed-mode and rotational remains similar to symmetric forcing, e = 0,with gradual increase in De threshold for a predominantly rotational re-sponse. For high forcing eccentricities, e > 0.5, all trajectories are in chaoticresponse. In general, energy ratio ER at a fixed e increases by decreasingDe.As we discussed in Chapter 2, we run experiments using four different60(𝑎𝑎) (𝑏𝑏)(𝑐𝑐)𝑓𝑓/𝑓𝑓𝑑𝑑 𝑓𝑓/𝑓𝑓𝑑𝑑𝑓𝑓/𝑓𝑓𝑑𝑑Figure 3.9: Example of Power spectrum for a. rotational, b. mixed-mode, and c. chaotic response. Red lines are from spectral anal-ysis of accelerometers data and black lines come from spectralanalysis of trajectory time series. The frequency is normalizedby the damped frequency of the column fd.columns with different elastic and damping properties, Table(2.4). Fig-ures 3.12, 3.13, 3.14, and 3.15 show the individual responses of columnsC1, C2, C3, and C4 to different forcing scenarios, respectively. Note, wehave defined a characteristics number for each column as γ = Tdωn, seeequation 2.6, representing the ratio of damping and natural oscillation timescales. In our experiments, the lowest value of γ = 8.15 belongs to columnsC1 and C4 while C1 has the shortest relaxation time, Td = 2.38 s, and C4has the longest oscillation time, 2pi/ωn = 3.85 s. The column C2 has theshortest oscillation time, 2pi/ωn = 1.44 s. The maximum value of γ = 14.161(𝑎𝑎) (𝑏𝑏)(𝑐𝑐)Figure 3.10: Power spectrum of all the experiments in a. rotational,b. mixed-mode, and c. chaotic response. The solid blackline represents the average value power spectrum for all theexperiments in each mode. The frequency is normalized bythe damped frequency of the column fd.and relaxation time, Td = 6.66 s, belongs to column C3.The response of the column with the shortest relaxation time (C1) toairflow with high eccentricities e > 0.5 is chaotic while in low eccentricitiese < 0.5 all three responses are observed (Fig. 3.12). Rotational solutionsare more favorite for high De and low eccentricity e while some mixed-modesolutions occur for mid-range values of eccentricity 0.15 < e < 0.45 andmid-range forcing intensity 7 × 103 < De < 1.8 × 104. For experiments at62Figure 3.11: Effects of forcing eccentricity e and intensity De on col-umn responses. Shapes represents responses as: © is rota-tional, ♦ is mixed-mode, and  is chaotic response. The col-orbar shows the value of ER.the area of e > 0.5 the energy ratio ER decreases by increasing De for aconstant e.Similar to column C1, all the responses of column with the shortestoscillation time (C2) for high airflow eccentricity e > 0.5 are chaotic andER decreases by increasing flow intensity De by fixing flow eccentricity e.Rotational solutions are observed in a limited forcing properties of 0.2 <e < 0.5 and De > 0.5 × 104. In contrast to column C1, all C2’s motion atsymmetry forcing e = 0 are chaotic.Rotational responses for the column with longest damping time (C3)occur only at e = 0, De = 1.4 × 104 and e = 0, De = 1 × 104 (Fig. 3.14)and mixed-mode responses are observed by slight increase in eccentricitye = 0.09, De = 1.3 × 104 or decrease in intensity of forcing e = 0, De =0.7× 104 compared respectively to rotational responses. ER increases fromits minimum at low eccentricity e = 0 and high intensity De = 1.4× 104 toits maximum at high eccentricity e > 0.7 and low intensity De < 0.3×104.63Figure 3.12: Effects of forcing eccentricity e and intensity De on C1responses. Shapes represents responses as: © is rotational, ♦is mixed-mode, and  is chaotic response. The colorbar showsthe value of ER.Figure 3.13: Effects of forcing eccentricity e and intensity De on C2responses. Shapes represents responoses as: © is rotational, ♦is mixed-mode, and  is chaotic response. The colorbar showsthe value of ER.64Figure 3.14: Effects of forcing eccentricity e and intensity De on C3responses. Shapes represents responses as: © is rotational, ♦is mixed-mode, and  is chaotic response. The colorbar showsthe value of ER.Similar to column C3, the column with longest oscillation time (C4) re-sponses chaotically to different airflows except at symmetry airflow e = 0 andhigh forcing intensity De = 1.6×104 which its response is rotational/mixed-mode and at e = 0.15 and De = 1× 104 which its response is mixed-mode.3.5 Discussion of experimental resultsIn this discussion we first overview different classes of responses and theirproperties. We then address: i. the significant role of airflow eccentricity ein dictating the class of column response; ii. the effect of De on the columnresponse for a symmetric airflow with e = 0; and iii. some unexpected resultsobserved in Fig. 3.11 where there are chaotic responses with symmetricairflow and high De.65Figure 3.15: Effects of forcing eccentricity e and intensity De on C4responses. Shapes represents responses as: © is rotational, ♦is mixed-mode, and  is chaotic response. The colorbar showsthe value of ER.3.5.1 OverviewIn this chapter, we first introduced rotational, chaotic, and mixed-moderesponses. Next we characterized different classes of responses based oncharacteristic differences in displacement, acceleration, and specific energytime series. A column in a rotational mode follows approximately constantorbits, has low energy ratios, ER < 0.1, and its energetic components arestatistically stationary. In contrast to rotational responses, a column witha mixed-mode response follows orbits that shift and change over a range oftime scales, and is marked by higher energy ratios than rotational mode,0.1 < ER < 0.2. The energetic content of top trajectory in a mixed-moderesponse vary over short and long time scales with some intermittent peri-ods of statistically stationary. A chaotic response mainly consists of lineardisplacements with large variations in the amplitude of oscillation. A widerange of energy ratios, ER > 0.2, is observed in a chaotic response whilecorresponding energetic components of top trajectory vary significantly over66time scales of ∼ Td.Vertical deflections along the column are distinct for each mode. Whereas,the deflection is monotonic in a rotational response (Fig. 3.1 b), it is spatiallyperiodic for a mixed-mode response (Fig. 3.3 b) and is complex during achaotic regime(Fig. 3.5 b). Comparing the time series of top displacement ofrotational, mixed-mode, and chaotic responses (Figs. 3.1 c,d, 3.3 c,d, 3.5 c,d)shows that a rotational response has the largest total amplitude of the de-flection, while the variability of the amplitude is minimized. On the otherhand, chaotic responses occur with the minimum amplitude of the deflectionand maximum variability of amplitude.Comparing spectral analyses of displacement time series at steady-stateregime (Figs. 3.9 and 3.10) with corresponding time series in damping regime(Fig. 2.7) shows that mechanical energy is transferred to the column at fre-quencies which are higher than the damped natural frequency of the column.To sustain the varied classes of motion in steady-state regimes, this energy istransferred to lower frequency, longer wavelength wagging (a reverse energycascade).We define the ratio Π of the power of the higher frequency second peak(f ≈ 10fd) to the first peak (f ≈ fd) corresponding to the damped naturalfrequency. Whereas for a rotational response Π ≈ 0.1, for mixed-mode andchaotic responses, Π ≈ 0.15 and Π > 0.3, respectively. This difference inΠ reflects the efficiency of energy exchange from higher frequency to lowerwagging frequency. That it is maximized for the rotational mode, consistentwith low energy ratio ER, indicating virtually all available kinetic energy istransferred to drive rotational orbits.3.5.2 Effects of forcing eccentricity e on column responseA key result from Figs. 3.11-3.15 is that wagging motions are always chaoticfor e > 0.5. Furthermore, for e < 0.5 , natural rotational modes are onlypreferred for High De. Thus, the key control over the character of waggingis e. We now discuss aspects of the underlying mechanics.Figure. 3.16 shows a column in symmetric airflow, if the column is in67resting position flow is distributed evenly adjacent to the column and theimposed restoring pressure force is equal in every direction. When the col-umn is displaced from its balance position to its left, the symmetric airflowis disturbed and because of the conservation of volume air speeds up at theleft side of the column while it decreases at the right side. Consequently,from conservation of energy expressed through Bernoulli equations appliedat specified heights on the left and right sidespl +12ρ0v2l = pr +12ρ0v2r , (3.1)where p and v represents flow pressure and velocity and subscripts l andr correspond to the left or right of the column, causes the pressure to de-crease at the left side and increase at the right side of the column. Thispressure difference introduces a net force to the left side (the Bernoulli ef-fect) and pushes the column further until elastic potential energy stored inthe column balances the flow-induced kinetic energy. An oscillation to theright occurs, in response, where the resulting elastic force drives the columnback beyond its initial resting position. Through this column movement,velocity increases at the right side and decreases at the left side, causing thenet pressure pushing the column to the right. Again, the column deflectsto the right until elastic potential energy in the column balances the flow-induced kinetic energy. The symmetric inflow implies a symmetric restoringpressure-force and leads to similar column displacements for both sides ofthe column.In this symmetric airflow condition if initial displacement of the columnoccurs in two perpendicular directions, two independent simultaneous os-cillations of the same frequency occur. The superposition of these simpleharmonic motions can produce elliptical or circular shapes depending onthe difference in amplitude and phase of the oscillations [14]. If the initialphase difference is pi/2 or 3pi/2, the motion is circular for the oscillationswith same amplitude. The initial phase difference of 0, pi, or 2pi results in alinear motion. For other values of phase differences an elliptical trajectoryis gained. Consequently, maintaining these oscillations produce rotational68responses.Where asymmetric airflow is imposed at the tank base, by contrast, theflow distribution around the column and resulting net pressure force at itsresting position is not uniform. In Fig. 3.17, the pressure at the right side ofthe column is higher than the column’s left side due to its lower flow velocitythrough equation 3.1. The resultant net pressure force pushes the columnto the left and the resting position of the column moves to the left. If thecolumn displaces from the new resting position, similar to symmetric airflowcondition, column will oscillate. Since the column oscillation center is on theleft of the channel centerline, the amplitude of column displacement fromits original resting position is greater in the left direction than the other.Extending asymmetric airflow condition problem to three dimensionscan result in different trajectories [14]. In this case, the superposition of topcolumn displacement at two perpendicular directions with amplitudes, andphases can lead to variety of trajectories.3.5.3 Effects of forcing intensity De on column responseWe have observed from time series and top trajectories (e.g., Fig.2.3) thatafter shutting down the airflow, the column is prone to oscillate in a ro-tational response. When the asymmetric airflow is exerted to the columnat the base, it introduces a linear net pressure force which can break therotational pattern leading to motions that include contributions from linearand rotational displacements. Consequently, column oscillates in a chaoticregime. By contrast, with a symmetric air flux forcing the net pressureforce is approximately zero (depending on physical noise). In this specialcase, neutralized radial pressure forces act to maintain rotational motions.However, we observe rotational, mixed-mode, and chaotic responses withsymmetric forcing e → 0. When e < 0.5, low airflow intensities (0.6 × 104)favor chaotic responses, and high De leads to an increasing probability ofmixed-mode (De < 104) and rotational modes (De > 104).To build understanding, we revisit the steady-state form of the equationof motion for a viscoelastic cantilever beam (equation 2.21) representing69RubbercolumnRubbercolumn𝑢𝑢𝑢𝑢 Gas velocityPressure variationRubbercolumnRubbercolumnInflowtankwallPlan view at baseFigure 3.16: Schematic of gas velocity and pressure variation aroundrubber column, and imparted pressure at the inlet for a sym-metric airflow condition. Longer arrows and darker colors rep-resent bigger magnitude of flow velocity and pressure forces.Note, imposed gas pressure force at base acting on a deflectedcolumn is always symmetric.70𝑢𝑢𝑢𝑢 + 𝑑𝑑𝑢𝑢RubbercolumnRubbercolumnRubbercolumntankwallInflowGas velocityPressure variationRubbercolumnPlan view at baseFigure 3.17: Schematic of gas velocity and pressure variation aroundrubber column, and imparted pressure at the inlet for an asym-metric airflow condition. Longer arrows and darker colors rep-resent bigger magnitude of flow velocity and pressure forces.Note, pressure variations imparted at base are not symmet-ric. Even if the column would prefer a rotational mode, thisinlet condition always perturbs the column away from axi-symmetry.71the experimental model. In the steady-state regime the equation of motionreduces to:∂4u(z, t)∂z4+ρAEI∂2u(z, t)∂t2=w′(z, t)EI(3.2)where w′ is the symmetric forcing driven as a result of difference betweenthe time-averaged and spatially uniform supply of mechanical energy to thecolumn and dissipation. We non-dimensionalize equation 3.2 by selecting thegap width δ as the length scale for deflection, column height L as the lengthscale for z direction along the column, T = L/V as the advective time scalefor gas flow where V is the average vertical velocity of the airflow. Definingu = u′δ, z = z′L, t = t′T , the equation of motion becomes (after droppingthe primes):∂4u∂z4+ Ω∂2u∂t2= wL4δ, (3.3)WhereΩ = (AL2I)(ρV 2E). (3.4)Ω is proportional to the ratio of total kinetic energy in the flow (∝ V 2) tothe elastic potential energy (∝ E), and to the Deborah number Ω ∝ De. Inthe limit of high De, equation 3.3 reduces to∂2u∂t2= wL4δΩ, (3.5)which represents column acceleration. In this case, the flow time scale isshorter than the time scale for the column to react, oscillations can growand a rotational response which corresponds to maximized amplitude of theoscillation occurs. In the limit of low De, by contrast, equation 3.3 reducesto∂4u∂z4= wL4δ, (3.6)which represents a bend in the column. In this case, the flow time scale islarger than the column reaction time scale and the kinetic energy deliveredto the column is stored as elastic potential energy that is not released overan advective time scale. Consequently, a chaotic response with e → 0 and72low De is a superposition of a series of arbitrary free oscillations of columnswith initial bending produced over T .3.5.4 Some unexpected results in e−De spaceFigure 3.13 shows that column C2’s response is chaotic when airflow is sym-metric e = 0 although we expect to have rotational responses for e = 0 andhigh intensity flows. To answer this contradiction, we consider that in orderto have symmetric airflow, not only the flow should be distributed evenlyat the base also the gap width between the column and hollow cylinder δwhen the flow is not introduced to the cylinder should be equal in all di-rections. Second investigations show that for most of the experiments runby column C2, the column was not standing straight in the center of thehollow cylinder, and δ was different for different directions. This also canexplain why rotational responses are observed at eccentricities in the rangesof 0.2 < e < 0.5; even though input flow is not perfectly symmetric, but theinitial deflection of the column cause flow to be distributed evenly aroundthe column and rotational responses are recovered.73Chapter 4Numerical SimulationsTo better understand and extend aspects of our experimental results relatedto the control of gas flux intensity and symmetry in forcing, we carry out2D Cartesian numerical simulations. In particular, consistent with previousstudies on the magma wagging model [6, 21] and the result in Chapter 3we observe two main frequencies: the natural frequency of column oscil-lation and the higher Bernoulli frequency related to the gas flux forcing.The steady-state regimes for rotational, mixed-mode, and chaotic modesimply also that energy is transferred from the Bernoulli mode to dampednatural frequency fd of wagging, confirming a hypothesis in [6]. Using two-dimensional numerical simulations to complement these results, we explorethe effects of total gas flux, physical column properties, and spatial asym-metry in the distribution of gas flux on the spectral properties of waggingbehavior.In this chapter, we first present the governing equations for a 2D elasticcolumn vibrating in response to an imposed gas flow and define related di-mensionless control parameters. We then briefly explain the numerical mod-eling strategy we implement using the COMSOL Multiphysics finite elementsoftware package. We present a case study to address the applicability of2D simulations for understanding wagging frequencies and driving Bernoullimode. Next, we investigate the effect of different defined dimensionlessparameters on the column deflection and pressure and velocity variations74around the column. Finally, we present a summary of the chapter.4.1 Numerical modelTo model our experiments, we implement the Fluid-Structure Interactionmodule of the COMSOL Multiphysics finite element software. The geom-etry in this study consists of two domains: fluid and solid mechanics do-mains (Fig. 4.1). The software provides continuity for the displacementand offers a conforming mesh at the coupled interface between flow andsolid domains. We solve the initial value-boundary value problem using theTime −Dependent solver to solve unsteady/time-transient problems. Theresults are post-processed with MATLAB.We restrict calculations to perfectly elastic column in 2D because thegoal is to capture the modes of oscillation generally, the inclusion of dampingprecludes modes with a period higher than damping time. We discuss thislimitation briefly at the end of the chapter.4.1.1 Geometry and parametersFigure 4.1 shows a schematic of the geometry of our numerical model. Anelastic column with a density of ρ, Young’s modulus of elasticity of E, andPoisson’s ratio of ν is fixed at the bottom and placed in the middle of achannel. When the column is at rest, the gap width between the columnand channel wall is constant and equal to δ at its right and left side. Airwith a density of ρ0 and viscosity of µ enters the channel with the averagevelocity of VR and VL from the right and left inlets and leave the channelfrom above where its pressure is balanced with the atmosphere. Table 4.1lists the parameters for each model.4.1.2 Equations and dimensionless parametersElastic column equations:As with the development of equation 3.2, if we neglect the viscous dampingterm in equation 2.21, to focus on the oscillation at steady-state regime, as75Extra fine FinerFine Normal𝛿𝛿𝐿𝐿𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑉𝑉𝐿𝐿𝑖𝑖𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃 = 0Elastic Column𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑉𝑉𝑅𝑅𝐷𝐷𝑥𝑥𝑧𝑧Figure 4.1: Schematic of the model used in COMSOL Multiphysics.An elastic column in dimensions of D×L is fixed at the bottom.The initial gap width between the elastic column and the sur-rounding wall is δ. Different types of mesh used in this study,with a number of elements, extra fine mesh: 93940, finer mesh:22550, fine: 15654, and normal: 9968.well as effects related to column buoyancy, the equation of motion for anelastic column is [4]:∂2u(z, t)∂t2+EIρA∂4u(z, t)∂z4= w(z, t), (4.1)where ρ is column’s density, A is cross-sectional area of column, E is Young’smodulus, I = piD4/64 is the second moment of inertia, u is deflection of thecolumn, z is the vertical distance along the column from the base, and76Parameter DescriptonL Column’s heightD Column’s diameterδ Gap widthρ Column’s densityE Column’s Elasticityρ0 Fluid densityµ0 Fluid viscosityVL Left inlet velocityVR Right inlet velocityTable 4.1: List of parameters set in the model.w is force per unit mass applied to the column. We non-dimensionalizeequation 4.1 by selecting the gap width δ as the length scale for x direction,column height L as the length scale for z direction, T = L/V as the advectivetime scale for gas flow where V = (VR + VL)/2 is the average of inputvelocities,. In contrast to our treatment of equation 3.2, we also introduceand explicit scale for pressure ρ0V2 that enables us to couple the responseof the column to the imposed gas flow we discuss next. Defining u = u′δ,z = z′L, t = t′T , w = w′ρ0V 2/(ρD), the equation of motion becomes (afterdropping the primes):∂2u∂t2+KG∂4u∂z4= w(z)1mε, (4.2)WhereKG = (EρV 2)(IAL2), (4.3)is column bending rigidity and indicates the importance of elastic stress (orenergy) to the dynamic pressure (kinetic energy) introduced from fluid flowto the column. In this development, note that KG = 1/Ω (equation 3.4).m = ρ/ρ0 is the density ratio between column and fluid, and  = δ/L andε = D/L are aspect ratios in our model. We explicitly non-dimensionalizedequation 4.1 in this way to enable m to emerge on the right hand side of the77equation as one of the functions of solid-fluid interaction.Flow equations:We solve the Navier-Stokes equation and continuity for the fluid phase inour model. The incompressible Navier-Stokes equation and the continuityequation can be written as:∂u¯∂t+ (u¯.∇)u¯− µ0ρ0∇2u¯ = − 1ρ0∇P, (4.4)∇.u¯ = 0, (4.5)where u¯ is the fluid velocity vector, P is the fluid pressure, ρ0 is the fluiddensity, µ0 is the dynamic viscosity, ∇ indicates the gradient differentialoperator, and ∇2 is the Laplacian operator. We non-dimensionalize equa-tion 4.4 and equation 4.5 by selecting V and V as the velocity scale in zand x direction, the gap width δ as the length scale for x direction, columnheight L as the length scale for z direction, T = L/V as the time scale, andρ0V2 for the pressure scale . In this case we define uz = u′zV where uz isflow velocity component in z direction, ux = u′xV where ux is flow velocitycomponent in x directions, x = x′δ, z = z′L, t = t′T , and P = P ′ρ0V 2.if we non-dimensionalize equation 4.4 and equation 4.5 and write them inseparate equations for z and x directions we have:∂uz∂t+ uz∂uz∂z+ ux∂uz∂x− Re(∂2uz∂z2+12∂2uz∂x2) = −∂P∂z, (4.6)2(∂ux∂t+ uz∂ux∂z+ ux∂ux∂x)− Re(2∂2ux∂z2+∂2ux∂x2) = −∂P∂x, (4.7)∂uz∂z+∂ux∂x= 0, (4.8)where Re = ρ0V δ/µ0 is the Reynolds number. If we assume  1 andRe 1, then equation 4.6 and equation 4.7 can be reduced to:∂uz∂t+ uz∂uz∂z+ ux∂uz∂x− 1Re(∂2uz∂x2) = −∂P∂z, (4.9)78Parameter equation C1 C2 C3 C4ε D/L 0.08 0.11 0.07 0.09 δ/L 0.006 0.008 0.011 0.014m ρ/ρ0 850 850 850 850Re ρ0V D/µ0KG EI/(ρAL2V 2)χ VR/VLTable 4.2: List of non-dimensional parameters set corresponding toexperiments.0 =∂P∂x. (4.10)Final set of equations needed to be solve by COMSOL to find elastic col-umn deflection and flow properties at each time are equations 4.2, 4.8, 4.9,and 4.10. On the basis of these equations, we use non-dimensional param-eters listed in Table 4.2 as inputs of each simulation. Also, we introduceχ = VR/VL as the ratio between two inlet velocities. Then COMSOL cal-culates the corresponding parameters set (Table 4.1) automatically usingnon-dimensional parameters. The value of each non-dimensional parametercorresponding to viscoelastic rubber column used in experiments is given inTable 4.2.4.1.3 Boundary conditionsWe define two Inlets at the bottom of the channel and introduce a parabolicvelocity profile with an average velocity of VR and VL for the right and leftinlets, respectively. The top of the channel is Outlet and is in atmosphericpressure (P0 = 0). The condition at the walls of the channel, as well as theelastic column, is no slip condition. The only limitation to elastic columnis a Fixed Constraint at its bottom which is attached to the channel (seeFig.4.1).79Figure 4.2: Top displacement of the column with 4 different meshsettings. In this model KG = 10, Re = 1000, m = 850, andχ = 0.5. The time required for each simulation was 45, 125,224, and 410 min for Normal, Fine, Finer, and Extra F inemesh settings, respectively.4.1.4 Mesh selectionWe use COMSOL’s built-in “physics mode” to optimize the mesh design.One of the most important factors to provide an accurate solution is thespatial resolutions of displacements within the column and the column-gasstress coupling. However, mesh density trades with the time required foreach simulation. Consequently, we ran a model with four different mesh sizesettings, Normal, Fine, Finer, and Extra F ine (see Fig. 4.2). Based onthe results, Finer settings and ExtraF ine settings results are similar, butsimulation with Finer setting is roughly two times faster than Extra F inesimulation. Our choice for the rest of the models will be Finer settings.4.1.5 Solvers and convergenceCOMSOL automatically chooses a solver suitable for the modeled physics.In our analysis, the solver is Fully Coupled with Direct Linear Solver. Thelinear solver was MUMPS, and all the internal parameters of the solver areset automatically.The time-step for time-dependent analysis was fixed as dt = T3/10 where80T3 = 2pi/ω3 is the period of third natural frequency of the elastic column.This time-step is short enough to capture all the oscillations occur with firstthree modes of oscillations. COMSOL will take additional steps in betweenthe defined time steps if necessary, e.g., CFL condition.4.2 ResultsIn this section, we first present the results of one specific simulation. Wethen do a parametric study to investigate how defined non-dimensional pa-rameters affect the model and solutions.For all the simulations we show in this study, we perturb the elasticcolumn from its original resting position in following way; we increase right(or left) inlet velocity from 0 at t = 0 to ≈ 2 × VR (or 2 × VL for the leftinlet) at t = 0.5×T , where T = L/V¯ and V¯ = (VR+VL)/2, then decrease itgradually to a steady value of VR (or VL) at t = T . We run each simulationfor 10T . Also in this set of simulations we only study dynamic dimensionlessvariables ofRe, χ, m, andKG, and aspect ratios remain constant as ε = 0.07,and  = 0.011.4.2.1 A case studyFigure 4.3 shows column displacement at different heights as a function oftime. The column vibrates with two distinct frequencies: a high-frequencyand low-frequency oscillations. The amplitude of both oscillations increasewith height above the column base. Whereas the high-frequency period is≈ 2T/7, the low-frequency oscillation period is ≈ 2 × T . In general, thecolumn tends to bend toward the left side where the average inlet velocityis two times greater than the average inlet velocity at the right side (notethat χ = VR/VL).Figure 4.4 shows the average velocity of the flow at two sides of theelastic column for locations with L/3, 2L/3, and L distant from the inlet.The same oscillations observed in column displacement are observed herewhile the oscillation amplitude is minimum at the bottom of the channel;it reaches its maximum value around the middle. It then decreases toward81𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝐿𝐿𝐿𝐿𝐿𝐿𝑅𝑅𝐿𝐿𝑙𝑙𝑙𝑙𝑙𝑙𝐿𝐿ℎ𝑖𝑖𝑖𝑖ℎFigure 4.3: Column displacement at L/3, L/2, 2L/3, and L distancesfrom the fixed side where L is length of the elastic column. Ar-rows and dash lines show the period of low- and high-frequencyoscillations. Non-dimensional parameters for this simulation areKG = 0.1, Re = 100, m = 850, and χ = 0.5. y axis is scaled bygap width δ and time is scaled by T = L/V .the top of the channel. Also, after inlet airflow becomes steady (t = T ), themaximum velocity in one side of the column coincides with the minimum ofthe velocity on the other side of the column.Figure 4.5 shows the pressure variation in the flow at two sides of theelastic column for locations with L/3, 2L/3, and L distant from the inlet.Pressure variations that are caused by flow velocity oscillations have thesame frequencies as velocity time series. The amplitude of oscillation inpressure time series decrease by moving from base of channel to column topon left or right sides. Pressure variations at a distance of L are minimalin amplitude, and it is very close to atmospheric pressure. Similar to thevelocity time series, the maximum of the pressure on one side of the columncoincides with the minimum of pressure on the other side of the column.82(𝑎)(𝑏)(𝑐)𝑓𝑙𝑜𝑤 𝑓ℎ𝑖𝑔ℎ𝑢𝑧(m/s)𝑢𝑧(m/s)𝑢𝑧(m/s)Figure 4.4: Average velocity of the airflow in right side and left sideof the elastic column at a. L/3, b. 2L/3, and c. L distancesmeasured from the bottom of the channel. Arrows and dashlines show the period of low and high-frequency oscillations.Parameters are KG = 0.1, Re = 100, m = 850, and χ = 0.5.Time is scaled by T = L/V .83(𝑎𝑎)(𝑐𝑐)𝑓𝑓𝑙𝑙𝑙𝑙𝑙𝑙 𝑓𝑓ℎ𝑖𝑖𝑖𝑖ℎ(𝑏𝑏)Figure 4.5: Pressure variations of the airflow in right and left sidesof the elastic column at a. L/3, b. 2L/3, and c. L distancesmeasured from the bottom of the channel. The solid and dashedblack lines present the pressure envelopes at the right and leftthe side of the channel, respectively. Arrows and straight dashlines show the period of low and high-frequency oscillations.844.2.2 Parametric studySimilar to forcing eccentricity e defined in Chapter 2, velocity ratio χ =VR/VL is a metric for the 2D asymmetry of the flow adjacent to both sidesof the elastic column. Reynold’s number Re is a metric for the availablekinetic energy transferred to the column by the imposed gas flow. TheBending rigidity KG is a metric for the importance of spring stiffness of thecolumn compared to the dynamic pressure of the flow, and density ratio mdefines the inertia of the column. In this section, we study the sensitivityof column response to these non-dimensional variables. To this end, in thefollowing figures, we change one of the variables and hold the rest fixed.Velocity ratio (χ):In this set of simulations non-dimensional parameters are as KG = 1,Re = 1000, m = 850, and we change velocity ratio between 0.3, 0.5, and 0.7.Note that χ = 1 is equivalent to symmetric forcing, and as χ → 0 repre-sents greater asymmetric forcing. Figure 4.6 a. shows column top displace-ment. There is only a low-frequency oscillation observed in displacementtime series. By decreasing χ, the amplitude of oscillation increases, but itsfrequency is not affected. Decreasing χ causes the column to bend furtherto the left, and the center of the oscillation moves further to the left side ofthe column.Figure 4.6 b. shows flow velocity variations and Fig 4.6 c. shows flowpressure variations at distance of 2L/3 from the bottom of the channel.Velocity time series show that by decreasing χ, the flow variation at twosides of the column increases, and the flow becomes more asymmetric.High frequency and low-frequency oscillation are observed at both ve-locity and pressure time series. Similar to high-frequency oscillation, theamplitude of low-frequency oscillation increases by decreasing the velocityratio χ while its frequency is intact.85𝑓𝑓𝑙𝑙𝑙𝑙𝑙𝑙𝑓𝑓ℎ𝑖𝑖𝑖𝑖ℎ(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝐿𝐿𝐿𝐿𝑓𝑓𝑅𝑅𝑢𝑢𝑧𝑧(m/s)Figure 4.6: a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the right chan-nel at 2L/3 from the bottom of the channel. Dashed lines inb. represent velocity at the left channel. Arrows and dashlines show the period of low and high-frequency oscillations.Non-dimensional parameters for these simulations are KG = 1,Re = 1000, m = 850, and χ changes for each simulation. Timeis scaled by T = L/V . 86Reynolds number Re:In this set of simulations non-dimensional parameters are as χ = 0.1, KG =1, m = 1000, and we change flow Reynold’s number Re between 10, 100,and 1000. Figure 4.7 a. shows column top displacement. We observe low-frequency oscillations in all three displacement time series, and the ampli-tude of the oscillation increases by increasing the Reynolds number. Also,the period of low-frequency oscillation slightly increases by increasing theRe number.Figure 4.7 b. shows flow velocity variations and Fig 4.7 d. shows flowpressure variations at distance 2L/3 from the bottom of the channel. Wemagnified velocity time series of cases with Re = 10 and Re = 100 by 100×and 10×, respectively, for visibility purposes. Also, pressure time series ofcases with Re = 10 and Re = 100 are magnified by 104× and 100×, respec-tively. The amplitude of oscillation for both high and low-frequency oscilla-tions in velocity and pressure time series increase by increasing Reynolds’snumber. The period of high-frequency oscillations does not vary by varyingReynolds number, while the period of low-frequency oscillation graduallyincreases by increasing Reynolds number.Bending rigidity KG:In this set of simulations non-dimensional parameters are as χ = 0.1, Re =1000, m = 1000, and we change bending rigidity KG between 1, 4, and10. Figure 4.8 a. shows column top displacement. Cases with smallerbending rigidity bend further, and the amplitude of oscillation increasesby decreasing bending rigidity KG. In the displacement time series, high-frequency oscillation is only observed when KG = 1, and the period oflow-frequency oscillations, broadly consistent with the results and analysisin Chapter 3. In addition to these spectral results, however, the amplitudeof the oscillations declines with increasing KG as well as their amplitudedecrease by increasing KG.Figure 4.8 b. shows flow velocity and Fig 4.8 d. shows flow pressurevariations at a distance of 2L/3 from the bottom of the channel. The am-87(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)𝑓𝑓𝑙𝑙𝑙𝑙𝑙𝑙𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝐿𝐿𝐿𝐿𝑓𝑓𝑅𝑅𝑓𝑓ℎ𝑖𝑖𝑖𝑖ℎ𝑓𝑓ℎ𝑖𝑖𝑖𝑖ℎ𝑢𝑢𝑧𝑧(m/s)Figure 4.7: a.Top displacement of the column, b. average flow ve-locity, and c. pressure variation of the flow in the right of thechannel at 2L/3 from the bottom of the channel. Parametersare KG = 1, χ = 0.1, m = 1000, and for different Re. Arrowsand dash lines show the period of low and high-frequency oscil-lations. Flow velocity is magnified 100× and 10× for Re = 10and Re = 100, respectively. Pressure variations for Re = 10and Re = 100 is 104× and 100× magnified, respectively.88plitude and period of both low and high-frequency oscillations decrease byincreasing bending rigidity. For the case with KG = 10, still, high-frequencyoscillations does not occur.To be consistent with the our high De experiments and to capture thehigh Ω limit identified in §3.5.3, We investigate the column response witha KG < 1 and ideally KG 1. However, in these limits, large columnaccelerations cause the column to interact with the walls over time scalesthat are sufficiently small to lead to the growth challenging numerical in-stabilities (violations of the CFL condition). One way to slow the growthrate of these instabilities and enable at least a qualitative exploration ofthis limit is to reduce Re to 1. Figure 4.9 shows a simulation where corre-sponding non-dimensional parameters are as χ = 0.9, Re = 1, m = 100, andKG = 0.01. In contrary to all previous simulations, the amplitude of theoscillation in column top displacement as well as in pressure and velocityvariations increase over time.Density ratio m:In this set of simulations non-dimensional parameters are as KG = 1,Re = 1000, χ = 0.5, and we change density ratio between 300, 600, and850. Figure 4.10 a. shows column top displacement. We observe that de-creasing m causes the column to bend further to the left, and the center ofthe oscillation moves further to the left side of the column (toward greateraverage flow velocity). The period of low-frequency oscillations increases bydecreasing the density ratio as simulation with m = 850 has the shortestperiod, and m = 300 has the longest period. In the displacement time series,high-frequency oscillations are much more significant when the density ratiodecreases.Figure 4.10 b. shows flow velocity, and Fig 4.10 d. shows flow pressurevariations at a distance of 2L/3 from the bottom of the channel. We observethat high-frequency oscillations are with the same period for all three cases,while by increasing density ratio, the amplitude of oscillation decreases.89(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝐿𝐿𝐿𝐿𝐿𝐿𝑅𝑅𝐿𝐿𝑙𝑙𝑙𝑙𝑙𝑙𝐿𝐿ℎ𝑖𝑖𝑖𝑖ℎ𝑢𝑢𝑧𝑧(m/s)Figure 4.8: a.Top displacement of the column, b. average flow veloc-ity, and c. pressure variation of the flow in the right channelat 2L/3 from the bottom of the channel. Arrows and dashlines show the period of low and high-frequency oscillations.Non-dimensional parameters for these simulations are χ = 0.1,Re = 1000, m = 1000, and KG changes for each simulation.Time is scaled by T = L/V .90(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)Figure 4.9: a.Top displacement of the column, b. average flow veloc-ity, and c. pressure variation of the flow in the right channel at2L/3 from the bottom of the channel. Non-dimensional param-eters for these simulations are χ = 0.9, Re = 1, m = 100, andKG = 0.01. Time is scaled by T = L/V .91(𝑎𝑎)(𝑏𝑏)(𝑐𝑐)𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝐿𝐿𝐿𝐿𝐿𝐿𝑅𝑅𝐿𝐿𝑙𝑙𝑙𝑙𝑙𝑙𝐿𝐿ℎ𝑖𝑖𝑖𝑖ℎ𝑢𝑢𝑧𝑧(m/s)Figure 4.10: a.Top displacement of the column, b. average flow veloc-ity, and c. pressure variation of the flow in the right channelat 2L/3 from the bottom of the channel. Arrows and dashlines show the period of low and high-frequency oscillations.Non-dimensional parameters for these simulations are KG = 1,Re = 1000, χ = 0.5, and m changes for each simulation. Timeis scaled by T = L/V .924.3 Discussions and summaryIn this discussion we first overview the response of an elastic column regard-ing to different dimensionless control parameters we defined at the beginningof this chapter. We then address: i. when do both high-frequency waggingand low-frequency Bernoulli mode appear?; ii. how can we relate the 2Dnumeric solution to the experiments in terms of different class of responses?4.3.1 OverviewIn this chapter, we studied the response of an elastic column to parallel flowin a channel with a Cartesian geometry as an approximation to our exper-iments. We first introduce four governing dimensionless parameters usingthe equation of motion for the elastic column and Navier-Stokes equations inthe fluid flow. Then, we investigate how these parameters affect the columnresponse using a case study and parametric studies.Through varied case studies, (Figs.4.3–4.5) we observe that the elasticcolumn oscillates with two distinct frequencies provided KG ≤ 1. In con-trast to the experiments, whereas high-frequency oscillations correspond tothe natural frequency of the column and low-frequency oscillation, whichcorresponds to Bernoulli’s effect of the flow adjacent to the column. Simula-tion results show that Bernoulli modes will always occur with the existenceof flow adjacent to the column. In contrast, high-frequency wagging modeis not always observed. Dimensional analysis of equation 4.2 shows that thecolumn’s period of oscillation grows with ∝ √1/KG. When KG increases,frequency of high-frequency oscillations increases, and higher amount of en-ergy is required to excite the oscillation. As a result, the mechanical energytransferred from low-frequency Bernoulli mode to high-frequency waggingmode has to be large enough to be able to excite wagging. Based on oursimulations, lowering the column’s bending rigidity (KG ↓), increasing theflow Reynolds number (Re ↑), increasing flow asymmetry (χ→ 0), and de-creasing the mass ratio between column and the flow (m ↓) can increase thechance of wagging modes to occur.In more details, we investigate Velocity ratio χ, Reynold’s number Re,93bending rigidity KG, and density ratio between column and gas m. Velocityratio χ = VR/VL is a metric that shows the asymmetric of the flow adjacentto the elastic column. Based on our definition, when χ = 1 flow is distributedevenly adjacent to the column, and the most asymmetric flow is when χ = 0.Varying velocity ratio χ does not change the frequency of high-frequency orlow-frequency oscillations but decreasing velocity ratio χ, enhancing the flowasymmetric, cause the amplitude of both oscillations to increase, see Fig.4.6.Investigation of Reynolds number Re, which is a metric for the availablekinetic energy in the flow to be transferred to the elastic column, shows thatperiod of high-frequency oscillations does not change by varying Reynoldsnumber while its amplitude increase if Reynolds number increases. On theother hand, not only the amplitude of low-frequency oscillation does increasewith increasing Reynolds number, but also the period of the low-frequencyoscillation increases with increasing Reynolds number, see Fig.4.7.Increasing density ratio m causes the period of low-frequency oscillationto decrease, but the period of high-frequency oscillations are constant. Andthe amplitude of both low and high-frequency oscillations increases whendensity ratio decreases, see Fig. 4.8.Bending rigidity KG represents the ratio of elastic stresses in the columnto dynamic pressure in the flow. Consistent with the discussion and analysisin §3.5.3, we find bending rigidity to be the only parameter that affects pe-riods of both low and high-frequency oscillation. Increasing bending rigidityreduces the period and amplitude of low and high-frequency oscillations inthe displacement time series as well as in velocity and pressure variations,see Fig. 4.8.Whereas the natural high-frequency oscillations of wagging are excitedfor KG ≤ 1 and are insensitive to all other parameters, the lower frequencyBernoulli mode depends on all the dimensionless parameters, and just itsvariation with the velocity ratio χ is negligible. Dimensional analysis ofequation 4.2 and equation 4.9 shows that the time scale of elastic columndisplacement and fluid flow variations is a function of Re, KG, and m.94Figure 4.11: Top displacement of the column. Non-dimensional pa-rameters for this simulation are KG = 0.01, Re = 1, χ = 0.9,and m = 100.4.3.2 The relation between the 2D numeric solution andexperiments in terms of different class of responses:As we discussed in Chapter 3, a rotational response can occur in a highDe limit or equivalently where Ω 1 (or KG 1). In this conditionthe column accelerates and oscillation continues. While in low De airflow(Ω 1) oscillation is a superposition of a series of arbitrary free oscillationand cause a chaotic solution.From the simulations, although the Bernoulli mode excites a high-frequencyoscillation where KG = 1, we observe that the only case that amplitude ofthe oscillation increases (column accelerates) over time is when KG = 0.01(Ω = 100). Figure 4.11 shows the top displacement of the column for asimulation with KG = 0.01, Re = 1, χ = 0.9, and m = 100 for a longer timethan in Fig. 4.9. In this case, the amplitude of oscillation increases for thewhole simulation. This case can potentially refer to a rotational responseas the flow is approximately symmetric (χ = 0.9), the amplitude of the col-umn displacement is nearly equal at both sides of the channel, and mostimportant, oscillation is growing not damping.95Chapter 5Implication for RealVolcanoesBercovici et al. [6] and Liao et al. [26] have shown that there is a correlationbetween column displacement inside a volcano and corresponding seismictremors. Assuming that the magma column inside the volcano oscillates withone of the rotational, mixed-mode, or chaotic responses, the related tremorwould have different seismic signatures. If there are a series of seismometersinstalled around a volcanic vent with the same radial distance, and themagma column is in a rotational mode, received tremor signals could becorrelated to each other with a constant lag [26]. If the magma column isin a mixed-mode response, seismic signals could be correlated with varyingtime lags in time. During a chaotic response, tremor will be observed withtransient periods of spatial correlation across the vent.Figure 5.2 shows how forcing exerted to a magma column inside an activevolcano can evolve from a constant and initially low eccentricity through arange of e−De conditions with pre-eruptive gas pressurization and throughthe eruption itself. Initially, Degassing through the annulus is plausiblyslow and uniformly distributed around the column (see Figs. 1.2 and 5.1)long times before an eruption happens.Gas flow forcing is at the low eccen-tricity indicated, and small Deborah number, De → 0. Consequently, themagma column is expect to oscillate in a chaotic mode with progressive gas-96Figure 5.1: Some snapshots of observed gas ring atthe edge of Karymsky volcano. Captured fromhttps://www.newsflare.com/video/227079/weather-nature/magnificent-ash-plumes-at-russias-karymsky-volcano97pressurization. Degassing becomes more intense over time. Increasing gasflow rate causes Deborah number to increase, De ↑, with the eccentricity un-changed. Through this evolution, column displacements transition smoothlyinto mixed mode. Transient spatial correlations in seismic signals of wag-ging across as well as around the volcanic event emerge and become morefrequent as an eruption approaches. As degassing rates increase towards apre-eruptive maximum, De becomes large and column displacements willenter an increasingly rotational regime. As a result, tremors become cor-related azimuthally, and a constant lag at different directions around thevolcanic vent can be observed.The smooth evolution at constant e changes with the first explosions thatdestroy parts of the annulus locally. The resulting spatial heterogeneity inpermeability causes asymmetry in the gas flux thus e to increase. In additionto natural wagging frequencies becoming spatially variable (see Chapter 1),the spatially complex gas flux forcing of these modes will drive wagginginto a chaotic mode. Consequently, the correlation between tremor signalsdecreases and finally disappears.98𝑒𝑒𝐷𝐷𝑒𝑒𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶 ≪ 𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝐶𝐶 ~ 𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐶𝐶𝐶𝐶𝐶𝐶𝐸𝐸Figure 5.2: Volcano’s activity in eccentricity e, Deborah De numberregime. A long time before the eruption, (t texplosion, vol-canic degassing is inefficient, De ≈ 0, and degassing is uniform,e ≈ 0. Closer to eruption time, degassing enhances, De ↑, whileannulus structure does not change. Eventually, enhancing de-gassing causes bubbles strain rate to increases above a criticalstrain rate, and fragmentation and explosion occur. Fragmenta-tion and explosion cause annulus to become heterogeneous andlead to non-uniform degassing around the volcano, e→ 1. Eachcolored area corresponds to the ideal condition for one of therotational, mixed-mode, or chaotic regimes based on the resultsin Chapter 3.99Chapter 6Concluding Remarks andFuture Research DirectionsThis thesis has investigated magma wagging model experimentally and nu-merically. The experimental setup and results are in Chapter 2 and Chap-ter 3, which studied Bernoulli mode and three distinct oscillation responses.Numerical modeling and results are presented in Chapter 4, which investi-gated Bernoulli and wagging modes in depth. These results are summarizedfirst in §6.1. Secondly, we discuss some limitations of the approach in thethesis, together with possible means of improvement (§6.2). The thesis closeswith suggestions for future research directions in this area (§6.3).6.1 SummaryIn this work, in Chapter 1 we have reintroduced magma wagging modelstudied in [6, 21, 26]. The main control parameters of the model, see equa-tion (1.6, are as following; the ratio of hydrostatic and dynamic gas pressureγ, the ratio of magma and gas densities β, ratio of the annulus thickness tothe magma column width λ, and the ratio of viscous forces arising in themagma column to the gas spring force in the annulus η.In Chapter 2 & Chapter 3, we have studied magma wagging model ex-perimentally to verify that the proposed mechanism of theory is a complete100description of the phenomena. Critically, to investigate whether the long pe-riod Bernoulli mode can continuously drive the predicted magma wagging.Also, experimental studies can help to identify new modes of oscillationincluding torsional and lateral wagging modes.We have introduced Deborah number De through equation 2.35 whichis a metric of airflow forcing in our experiments. Also, Deborah number isrelated to rate of available energy in the airflow to the rate of dissipation ofenergy in the column. It is analogues to η and captures the basic mechanicalproperties of magma wagging model.The other important feature that the experiments investigate is the dis-tribution of non-uniform gas flux around the column. The non-uniform gasflux is related to explosion and fragmentation. The fragmentation mod-ifies the annulus structure and changes permeability of the annulus. Weextend the model by introducing heterogeneous gas flux around the column.We have defined forcing eccentricity e, equation 2.34, as the parameter forcharacterizing the asymmetry of the airflow around the column.We have used energy components of the trajectory of top of the columnto characterize the response of the elastic column to airflow forcing by usingdimensionless energy ratio ER equation 2.49. Then, we have mapped theexperimental results on a De−e parameter space and investigate the favoriteconditions for each of the rotational, mixed-mode, and chaotic responses.In Chapter 4, we have studied magma wagging numerically in a 2DCartesian geometry. We have derived the control parameters of Re, m,and KG, through governing equations of solid and fluid domains. We haveextended the study by adding velocity ratio χ between average flow speedat two sides of the column as a metric similar to forcing eccentricity e tothe problem. We then investigate the column response to variety of controlparameters in terms of frequency and amplitude of the oscillation.Here are the main conclusion from our study:Experimental results:• Displacements in the unforced damping regime confirm essential spec-101tral predictions made by [21] and [6] that wagging oscillation occurswith an amplitude that declines over time scale of Td. These exper-iments also show that the preferred mode of unforced oscillation is arotational displacement.• During steady-state regime, there is a relatively high frequency Bernoullimode which occurs along with lower frequency wagging. Mechanicalenergy is transferred from high frequency Bernoulli mode to lower fre-quency wagging to sustain oscillations through a reverse energy cas-cade.• Rotational responses in a symmetric airflow forcing e = 0 with highDe, and in damping regime, confirms qualitatively the model presentedby Liao et al.[26] that column goes under different whirling motions.• Rotational modes only occur with symmetric airflow e→ 0 for stronggas forcing, high De, where flow timescale is shorter than columnresponse time (Ω 1) and column accelerates.• We extend the predicted motion by [26] for a symmetric airflow e = 0and find new classes of chaotic and mixed-mode responses correspond-ing to low and intermediate De, respectively.• Chaotic responses with e = 0 and low De are a consequence of i. lackof total kinetic energy to stabilize the massive column into rotational,or ii. larger flow timescale than column reaction time scale leads to aseries of arbitrary free oscillation. The superposition of these arbitraryfree oscillations will be a chaotic response.• When airflow eccentricity increases as e → 1, trajectories becomeincreasingly chaotic, and rotational responses are impossible wheree > 0.5.• Chaotic responses for e > 0.5 are driven by lateral strong pressuregradients or non-zero net pressure forces related to these boundaryconditions. Non-zero net pressure forces or lateral pressure gradients102drive linear motions that perturb displacements that are otherwiserotational.Numerical results:• A relatively low frequency Bernoulli mode is observed that providesthe energy for higher frequency wagging mode (energy cascade fromlower frequency to higher frequency oscillations, similar to predictionsin Bercovici et al.[6]).• The frequency of the wagging mode is only a function of column prop-erties (KG) and wagging frequency increases by increasing KG.• The amplitude of the wagging mode depends on coupled physics offluid flow and solid mechanics. By lowering the column’s bendingrigidity (KG ↓), increasing the flow Reynolds number (Re ↑), increas-ing flow asymmetry (χ → 0), and decreasing the mass ratio betweenthe column and the flow (m ↓) the amplitude of wagging increases.• Low frequency Bernoulli mode is a function of coupled physics betweenfluid flow and the elastic column. Increasing Re, decreasing KG, anddecreasing m cause the the frequency of Bernoulli mode to increase.• Simulations with very small KG (equivalent to very large Ω) confirmsthat the acceleration can occur in column displacement and probablylead to a rotational solution.6.2 Limitations of the experimentsAnalog experiments of magma wagging are challenged by scaling issues wediscuss in Chapter 2. Although we have made several advances doing theexperiments and numerical modelings, we must also acknowledge some im-provements that can be added:• The permeable vesicular annulus can vary in space, particularly duringan eruption. We modeled this heterogeneity with imposing asymmetricairflow, however, the elasticity of the column is assumed homogeneous.103This assumption is not necessarily correct and non-homogeneous elas-tic property can be included in the model. It can affect the springforcing which results in different wagging frequency.• A potential weakness is that the column is viscoelastic in our exper-iments and more likely viscoplastic in real life with elasticity on theoutside. However, the magma column can be modeled as a highly vis-cous fluid encircled by an elastic layer. This design leads to anothercontrol parameter, viscosity of the fluid, which can be easily adjusted.• The magma column buoyancy forces are not included in the originalmodel, while it has influenced our experiments as we observed someunexpected behaviors, see more details in §3.5.4. Similar to our exper-iments that the column wants to bend because of the gravity, magmacolumn wants to spread out into the annulus. This problem can beinvestigated in more details.6.3 Future research directions: exploring volcanoseismic dataThis study can be extended to a real volcano by testing the predictions madein Chapter 5. If there is a well instrumented volcano to monitor tremor ac-tivity, received tremor signals can be studied and correlation between themcan be monitored. For example, there are a series of seismic stations aroundthe crater of Redoubt volcano. And the tremor signals related to recentactivities are available from these station. Bercovici et al.[6] studied thetremor signals measured at two station across the Redoubt vent prior to its2009 eruption and showed that the signals are out of phase by approximatelyhalf the dominant oscillation period. Later, Liao et al.[26] studied the crosscorrelation between the tremor signals from pairs of stations around Re-doubt vent and suggested a circular wagging that can explain the observedradiation pattern.As an extension to these studies, the evolution of magma column motionat early stages of volcanic activity, prior to the eruption, and during the104eruption can be investigated through studying the correlation of tremorsignals around an active volcano. If tremor signal are correlated betweenstations with a constant time lag, the magma column is wagging in rotationalmode. Columns wagging in the mixed-mode response will be characterizedby transient correlated tremor among stations with time-varying time lags.Chaotic wagging is characterized by no correlation of tremor signals amongdifferent stations.105Bibliography[1] K. Aki and R. Koyanagi. Deep volcanic tremor and magma ascentmechanism under Kilauea, Hawaii. J. Geophys. Res., 86(B8):7095,1981. → page 2[2] N. S. Bagdassarov, D. B. Dingwell, and S. L. Webb. Viscoelasticity ofcrystal- and bubble-bearing rhyolite melts. Phys. Earth Planet. Inter.,83(2):83–99, 1994. → page 7[3] C. K. Batchelor and G. K. Batchelor. An Introduction to FluidDynamics. Cambridge University Press, 2000. → page 45[4] C. F. Beards. Engineering Vibration Analysis with Application toControl Systems. Elsevier, 1995. → page 76[5] J. P. Benoit and S. R. McNutt. New constraints on source processes ofvolcanic tremor at arenal volcano, costa rica, using broadband seismicdata. Geophys. Res. Lett., 24(7194):449–452, 1997. → pages 2, 15[6] D. Bercovici, A. M. Jellinek, C. Michaut, and D. C. Roman. Volcanictremors and magma wagging: gas flux interactions and forcingmechanism. Geophys. J. Int., 195(2):1001–1022, 2013. → pagesiii, xiii, xiv, 5, 7, 10, 11, 12, 13, 14, 15, 16, 17, 18, 22, 44, 74, 96, 100, 102, 103, 104[7] B. A. Chouet. Long-period volcano seismicity: its source and use ineruption forecasting. Nature, 380(6572):309–316, 1996. → pages 2, 4[8] L. Collier, J. W. Neuberg, N. Lensky, V. Lyakhovsky, and O. Navon.Attenuation in gas-charged magma. J. Volcanol. Geotherm. Res., 153(1-2):21–36, 2006. → pages xii, 6[9] R. S. Crosson and D. A. Bame. A spherical source model for lowfrequency volcanic earthquakes. J. Geophys. Res., 90(B12):10237,1985. → page 4106[10] R. P. Denlinger and R. P. Hoblitt. Cyclic eruptive behavior of silicicvolcanoes. Geology, 27(5):459–462, 1999. → page 4[11] C. Donaldson, C. Caudron, R. G. Green, W. A. Thelen, and R. S.White. Relative seismic velocity variations correlate with deformationat Kı¯lauea volcano. Sci. Adv. , 3(6):e1700219, 2017. → page 9[12] D. Dzurisin. A comprehensive approach to monitoring volcanodeformation as a window on the eruption cycle. Rev. Geophys., 41(1),2003. → pages 1, 2[13] J. C. Eichelberger, C. R. Carrigan, H. R. Westrich, and R. H. Price.Non-explosive silicic volcanism. Nature, 323(6089):598–602, 1986. →page 5[14] A. P. French. Vibrations and Waves. CRC Press, 2017. → pages68, 69[15] H. M. Gonnermann and M. Manga. The fluid mechanics inside avolcano. Annu. Rev. Fluid Mech., 39(1):321–356, 2007. → pagesxii, 6, 9[16] A. Goto. A new model for volcanic earthquake at Unzen Volcano:Melt Rupture Model. Geophys. Re. Lett., 26(16):2541–2544, 1999. →page 4[17] M. T. Gudmundsson, R. Pedersen, K. Vogfjo¨rd, B. Thorbjarnardo´ttir,S. Jakobsdo´ttir, and M. J. Roberts. Eruptions of Eyjafjallajo¨kullVolcano, Iceland. EOS, Trans. American Geophys. Union, 91(21):190,2010. → page 1[18] M. Hellweg. Physical models for the source of Lascar’s harmonictremor. J. Volcanol. Geotherm. Res., 101(1):183 – 198, 2000. → page4[19] R. E. Holasek, S. Self, and A. W. Woods. Satellite observations andinterpretation of the 1991 Mount Pinatubo eruption plumes.J. Geophys. Res.: Solid Earth, 101(B12):27635–27655, 1996. → page 1[20] A. J. Hotovec, S. G. Prejean, J. E. Vidale, and J. Gomberg. Stronglygliding harmonic tremor during the 2009 eruption of RedoubtVolcano. J. Volcanol. Geotherm. Res., 259:89–99, 2013. → pagesxi, 2, 3107[21] A. M. Jellinek and D. Bercovici. Seismic tremors and magma waggingduring explosive volcanism. Nature, 470(7335):522–525, 2011. →pages iii, xii, 4, 5, 6, 7, 9, 10, 17, 22, 44, 74, 100, 102[22] J. B. Johnson and J. M. Lees. Plugs and chugs-seismic and acousticobservations of degassing explosions at Karymsky, Russia and Sangay,ecuador. J. Volcanol. Geotherm. Res., 101(1-2):67–82, 2000. → pagesxiii, 4, 15, 16[23] C. Klug and K. V. Cashman. Permeability development invesiculating magmas: implications for fragmentation. Bull. Volcanol.,58(2-3):87–100, 1996. → pages xii, 6[24] K. I. Konstantinou and V. Schlindwein. Nature, wavefield propertiesand source mechanism of volcanic tremor: a review.J. Volcanol. Geotherm. Res., 119(1-4):161–187, 2003. → page 2[25] R. C. Leet. Saturated and subcooled hydrothermal boiling ingroundwater flow channels as a source of harmonic tremor.J. Geophys. Res.: Solid Earth, 93(B5):4835–4849, 1988. → page 4[26] Y. Liao and D. Bercovici. Magma wagging and whirling: excitation bygas flux. Geophys. J. Int., 215(1):713–735, 2018. → pagesiii, xiv, xv, 17, 19, 20, 21, 22, 45, 96, 100, 102, 104[27] Z. Lu and D. Dzurisin. Ground surface deformation patterns, magmasupply, and magma storage at okmok volcano, alaska, from InSARanalysis: 2. coeruptive deflation, july–august 2008. J. Geophys. Res.,115, 2010. → pages 2, 9[28] S. R. McNutt. Volcanic seismology. Annu. Rev. Earth Planet. Sci., 33(1):461–491, 2005. → page 2[29] S. R. McNutt and T. Nishimura. Volcanic tremor during eruptions:Temporal characteristics, scaling and constraints on conduit size andprocesses. J. Volcanol. Geotherm. Res., 178(1):10–18, 2008. → page 2[30] P. A. Nadeau, J. L. Palma, and G. P. Waite. Linking volcanic tremor,degassing, and eruption dynamics via SO2 imaging.Geophys. Res. Lett., 38(1):L01304, 2011. → pages xiv, 15, 17[31] J. Neuberg. Characteristics and causes of shallow seismicity inandesite volcanoes. Phil. Trans. Math. Phys. Engin. Sci., 358(1770):1533–1546, 2000. → page 2108[32] S. Okumura, M. Nakamura, and A. Tsuchiyama. Shear-inducedbubble coalescence in rhyolitic melts with low vesicularity.Geophys. Res. Lett., 33(20), 2006. → page 5[33] S. Okumura, M. Nakamura, A. Tsuchiyama, T. Nakano, andK. Uesugi. Evolution of bubble microstructure in sheared rhyolite:Formation of a channel-like bubble network. J. Geophys. Res., 113(B7), 2008. → page 5[34] A. C. Rust, N. J. Balmforth, and S. Mandre. The feasibility ofgenerating low-frequency volcano seismicity by flow through adeformable channel. Geol. Soc. Lond. Spec. Publ., 307(1):45–56, 2008.→ page 4[35] R. S. J. Sparks. Dynamics of magma degassing.Geol. Soc. Lond. Spec. Publ., 213(1):5–22, 2003. → pages 2, 5[36] H. Tuffen, R. Smith, and P. R. Sammonds. Evidence for seismogenicfracture of silicic magma. Nature, 453:511–514, 2008. → page 9[37] K. Unglert and A. M. Jellinek. Feasibility study of spectral patternrecognition reveals distinct classes of volcanic tremor.J. Volcanol. Geotherm. Res., 336:219–244, 2017. → page 2[38] S. R. Young, R. S. J. Sparks, W. P. Aspinall, L. L. Lynch, A. D.Miller, R. E. A. Robertson, and J. B. Shepherd. Overview of theeruption of Soufriere Hills Volcano, Montserrat, 18 July 1995 toDecember 1997. Geophys. Res. Lett., 25(18):3389–3392, 1998. → page1109
Thesis/Dissertation
eng
Geophysics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
An experimental study of volcanic tremor driven by magma wagging
Text