Michael Reinwald
email: michael.reinwald@upmc.fr

Quentin Grimal

Stefan Catheline

Lapo Boschi

Super-resolution in near-field acoustic

teaching and research institutions in France or abroad, or from public or private research centers.

1 Introduction

26

It is well known [START_REF] Arthur H Benade | Fundamentals of musical acoustics[END_REF] that anatomy contributes to the 27 task of auditory source localization, as its effects on an 28 acoustic signal, described by the head-related transfer 29 function (HRTF) [START_REF] Wenzel | Localization using nonindividualized head-related transfer functions[END_REF][START_REF] Hu | A local representation of the head-related transfer function[END_REF], can be seen as a spectral fil-30 ter and depend on the location of the signal's source.

Human auditory source localization mostly relies on 32 differences in the phase and amplitude of signals per-33 ceived by the two ears, as well as "spectral cues", or frequency-dependent effects associated with the shape 35 of the pinnae and, possibly, other features of the body 36 ([34]).

37

Building on the work of Catheline et al. [7], we ex-plore here the specific role of elastic waves mediated in a skull-shaped object mimicking bone-conducted sound. While this study does not address the issue of whether and how bone conducted sound is employed by the human auditory (ears/brain) system, our goal is to determine whether these reverberated signals contain specific information about the reconstruction of the position of an auditory source, especially in the near-field. This could be relevant to current efforts in the study of bone conduction sound [START_REF] Jack J Wazen | Results of the bone-anchored hearing aid in unilateral hearing loss[END_REF][START_REF] Bruce N Walker | Evaluation of bone-conduction headsets for use in multitalker communication environments[END_REF][START_REF] Stenfelt | Implantable Bone Conduction Hearing Aids[END_REF][START_REF] Littler | Hearing by bone conduction and the use of bone-conduction hearing aids[END_REF][START_REF] Stenfelt | Boneconducted sound: physiological and clinical aspects[END_REF]. Using the principle of acoustic time reversal [START_REF] Fink | Time-reversed acoustics[END_REF][START_REF] Fink | Time-reversal acoustics in complex environments[END_REF], we convert the signal recorded by two receivers into the spatial coordinates of a source in the horizontal and sagittal plane, and evaluate the resolution with which the source position is thus reconstructed. Catheline et al. [7] showed via a time reversal experiment with a dry skull that in-skull elastic wave propagation provides information about spatial positioning of a sound source. They found that their time reversal algorithm, using elastic waves alone, received at two or only one recording transducer mimicking the ear, successfully reconstructed the source position(s), for single as well as multiple sources. The spatial resolution of this method was found to decrease with increasing distance between the skull and the sound source. This is in good agreement with the far-field diffraction law, which provides a relationship between the spatial resolution and the distance separating the antenna (skull) from the source. Our objective is to expand the early work of Catheline et al. [7] and Ing et al. [START_REF] Ing | In solid localization of finger impacts using acoustic time-reversal process[END_REF] to [START_REF] Arthur H Benade | Fundamentals of musical acoustics[END_REF] analyze the resolution of the same algorithm for a skull-shaped antenna specifically in the near-field, i.e., the sound source is placed closer than one wavelength to the skull, and (2) to evaluate the directionality of the algorithm, i.e. evaluate changes in resolution with respect to angular position of the sound source.

In this study, we conduct a suite of experiments on a simple setup, equivalent to the setup used in Catheline et al. [7], consisting of two recording transducers glued to a replica human skull. Sound is generated by a small speaker deployed at a variety of distances and azimuths. Our results show in particular that, in the near-field, the resolution with which we reconstruct the source position changes as a function of azimuth with respect to the skull and is clearly influenced by complex features of the skull such as the orbital cavities. Furthermore we achieve superresolution throughout all angles for sources very close to the skull.

Similarly minded experiment have been conducted

in recent years e.g. in the context of optics, where imaging with evanescent waves allows to surpass the classical diffraction limit; the super-resolution of nearfield microscopes is piloted by their probe size [START_REF] Dieter W Pohl | Optical stethoscopy: Image recording with resolution λ/20[END_REF][START_REF] Lewis | Development of a 500 å spatial resolution light microscope: I. light is efficiently transmitted through λ/16 diameter apertures[END_REF].

In this context, a source [START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy[END_REF][START_REF] Robert M Dickson | On/off blinking and switching behaviour of single molecules of green fluorescent protein[END_REF]4] or scatterers [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep superresolution vascular imaging[END_REF] smaller than a wavelength, placed within the medium can be detected in the far-field with super-resolution as well. Time reversal experiments can also surpass the diffraction limit when resonators are placed near a source [START_REF] Lerosey | Focusing beyond the diffraction limit with far-field time reversal[END_REF][START_REF] Rupin | Super-resolution experiments on lamb waves using a single emitter[END_REF] or when an acoustic sink is used [9]. To a lesser degree, near-field details can sometimes be extracted from the far-field using sophisticated algorithms such as inverse filter [8] or MU-SIC [START_REF] Simonetti | Localization of pointlike scatterers in solids with subwavelength resolution[END_REF]. Experiments with metamaterials, superlenses and hyper-lenses [START_REF] Brian | Negative refraction makes a perfect lens[END_REF] demonstrate moderate sub-diffraction imaging down to a quarter of the optical wavelength. All these techniques use different terminology but they all require some near-field measurements.

Because very few studies in psychoacoustics have explored human sound localization performances for nearby sources [START_REF] Parseihian | Reaching nearby sources: comparison between real and virtual sound and visual targets[END_REF], we are unable to determine whether the resolution achieved by our algorithm reproduces the performance of human listeners using bone conducted sound. While we do find that elastic waves contain sufficient information to successfully reconstruct source positions in the near-field, we cannot yet establish whether a similar capability is achieved by the human auditory system.

Methods

The experimental setup is based on the previously conducted experiment of Catheline et al.

149

A sketch of the experimental setup in the horizontal 150 plane is shown in Figure 1.

151

We checked that the sensors solely measure the vi-152 bration of the skull and are unresponsive to airborne 153 sound. This ensures that the time reversal algorithm 154 will utilize only elastic waves. Additionally, the in-155 fluence of the foam platform used to place the loud-156 speaker at certain distances has been tested to have 157 no influence on sound emission of the loudspeaker.

158

The first part of the experiment consists of recording the signals at the sensors for each speaker position. The speaker emits a chirp signal c(t) with a duration of 1 s and a linear frequency distribution between 0 Hz and 6 kHz. The function in time for such a chirp of duration T , minimum frequency f 0 and maximum frequency f 1 reads

c(t) = sin Φ 0 + 2π f 0 t + k 2 t 2 , (1)
with the initial phase Φ 0 at time t = 0 and the chirpy-159

ness k = f1-f0
T (in our case k = 6000Hz/s), also 160 known as the rate of frequency range across the chirp.

161

For each distance D the source positions in the hori-162 zontal plane are defined by the azimuth ϕ.

The recorded signal s at one of the sensors' location r, writes

s(ϕ 0 , r, t) = c(t) * G(ϕ 0 , r, t), (2)
where * denotes convolution, ϕ 0 is the source position Following Fink [START_REF] Fink | Acoustic Time-Reversal Mirrors[END_REF], the received signal s(ϕ 0 , r, -t) is time-reversed, i.e flipped with respect to time. It must then be backward propagated to any possible location ϕ i . This is equivalent to convolving s(ϕ 0 , r, -t) with the Green's function G(ϕ i , r, t). Since we do not have access to G(ϕ i , r, t), but we do have a library of recordings of s(ϕ i , r, t) for all possible values of ϕ i , we implement

T i (ϕ 0 , r, t) = s(ϕ 0 , r, -t) * s(ϕ i , r, t) = = c(-t) * G(ϕ 0 , r, -t) * c(t) * G(ϕ i , r, t). (3)
The term G(ϕ 0 , r, t) * G(ϕ i , r, -t) is the transfer function of such a time reversal algorithm and, in terms 187 of signal analysis, represents a matched filter [START_REF] Fink | Acoustic Time-Reversal Mirrors[END_REF].

188

This convolution coincides with the cross-correlation of G(ϕ 0 , r, t) and G(ϕ i , r, t) ( [START_REF] Draeger | One-channel time-reversal in chaotic cavities: Theoretical limits[END_REF][START_REF] Derode | Recovering the green's function from field-field correlations in an open scattering medium (l)[END_REF]). For each source position ϕ 0 , the signal processing procedure consists of implementing Equation 3, i.e. analytically crosscorrelating the signals, and of finding the maximum value, with respect to time, of the time-reversed wave field T i for each ϕ i . The resulting function F (ϕ i ) is dubbed "spatial focusing function" (shortly, focusing function), as this procedure is equivalent to evaluating whether (and with what resolution) the timereversed and backward-propagated wave field is able to reconstruct the original source position ϕ 0 . The focusing function is next normalized with respect to its maximum; It is then reasonable to assume that, the closer F (ϕ i ) is to 1 (i.e., identical Green's functions) for a given value of ϕ i , the closer ϕ i is to the original source ϕ 0 . This method can be interpreted as a pattern recognition system, that identifies, from an acoustic reference library, the Green's function corresponding to the actual position of the source, and so determines the position of the source.

The invariance under time reversal is lost if the propagation medium has frequency-dependent attenuation. This introduces a first-order time derivative in the governing propagation equation. However, the theorem of spatial reciprocity is still valid, i.e. there is a loss of amplitude in the time-reversed vs. forward propagating wave field, but this does not affect source-localization resolution (does not affect the location of the focus of the time-reversed wave field) provided that signal-to-noise ratio of recorded data is sufficiently high. We have accordingly chosen to carry out our experiments at frequencies that are well caught by our receiving system.

We take both sensors into account by computing the mean of the focusing functions of the two signals.

In order to investigate the role of different frequency contents, the originally measured signals are successively filtered with varying low-pass filters with maximum frequency f max . Following e.g. [START_REF] Ing | In solid localization of finger impacts using acoustic time-reversal process[END_REF]5,[START_REF] Tsogka | Time reversal through a solid-liquid interface and superresolution[END_REF], we estimate the spatial resolution of our time reversal algorithm by analyzing the -3 dB width ∆p of F (ϕ i ) for each given source position (angle ϕ and distance D between the source position and the skull) and various smallest wavelengths λ min = c/f max (with c = speed of sound in air).

We compare our resolution estimates against the apparent aperture A of our skull-shaped antenna, as defined by Catheline et al. [7], through the far-field diffraction law

A = D • λ min 2∆p . (4)
While resolution as defined here is known to follow the diffraction-law in the far-field [7], that is not the case in the near-field, where Equation 4 is only used here for the sake of comparison. diameters of the skull, close to 10 and 15 cm, in the horizontal and sagittal planes, respectively.

The measurement points in the near-field (at distances smaller than one wavelength) lie on the same linear fit (i.e. same apparent aperture) as the points for measurements in the far-field although Equation 4does not hold true in the near-field. In the near-field, i.e. for sources closer than one minimum wavelength away from the skull, source positions can still be resolved with the same angular resolution which results in super-resolution in space, i.e. -3 dB widths below 0.5 λ min (see Figure 4). While one could infer that the diffraction limit also holds true in the near-field, our results are purely empirical; any values below the previously formulated diffraction limit are not represented in Equation 4. We speculate that they can be ascribed to the near-field contribution of evanescent waves.

Our far-field data is in agreement with Equation 4 and the previous findings of Catheline et al. [7]. In addition, we are able to achieve the same angular resolution as stated in the far-field diffraction law in the near-field (sound sources at below-wavelength distances) leading to super-resolution.

Directional variation in resolution

We furthermore investigate the directional variation of resolution of the time reversal analysis in the horizontal plane. The angular variations in resolution of our time reversal scheme in the near-field are visualized in Figure 5 showing the values of A (top) and ∆p (bottom) with respect to the source azimuth ϕ for different source distances (5 cm, 12 cm and 20 cm and 100 cm). All data is filtered to have a maximum frequency of 3 kHz. The reason for an offset of around 2 -3 • to the center (ϕ = 0 • ) is due to a limited accuracy in the manual placement of the center position and the center of the rotation axis.

In the far-field, the apparent aperture does not vary with azimuth (see 100 cm data in Figure 5) and is ponentially over roughly one wavelength [START_REF] De | Focusing prop-445 erties of near-field time reversal[END_REF] and can 338 be created at a boundary between two media through certain incident angles of a propagating wave [START_REF] Fink | Time reversal of ultrasonic fields. i. basic principles[END_REF]. Usually, their effect is not measured in the far field and the far-field diffraction law (Equation 4) does not account for such effects, limiting the resolution of time reversal. However, if near-field components of the wavefield are measured and incorporated in the time-reversal algorithm, subwavelength information, that is carried by evanescent waves, is incorporated in the time-reversal process, leading to super resolution [START_REF] Lerosey | Focusing beyond the diffraction limit with far-field time reversal[END_REF].

All these results are also approximately achieved via a one-sided evaluation of the signals, i.e. when only one receiver is used.

In summary, our data shows large variations in resolution in the near-field, depending on the position of the source relative to the geometric complexities of the skull.

Conclusion

In this study we measured elastic wave signals in a replica of a human skull due to an incident airborne sound emitted by a source at various distances and orientation with respect to the skull. Our goal was to investigate the physical limits of a sound-localization algorithm that uses full waveform information and the information contained in elastic waves propagating in the skull bone. While we do not at all claim to directly reproduce the sound localization "algorithm" that exists in the human ear-brain system, our quantification of these limits may be considered as a point of comparison in near-field psychoacoustics experiments.

We showed that the resolution of a time reversal scheme using a skull-shaped antenna with one or two receivers is consistent with the diffraction law in the far-field. The apparent apertures in the horizontal and sagittal planes are roughly consistent with the horizontal and vertical extent of the skull. Interestingly, the apparent aperture in the near-field is markedly increased (more than 3 times its value in the far-field) in the horizontal plane and at specific angles. In that case we can achieve super-resolution that may be associated to the non-negligible contribution of evanescent waves in the near-field.

Our results suggest that anatomical details of the skull give rise to complex features of the radiated sound field in the near-field, enabling sub-wavelength focusing and directional changes in resolution. We clearly find the influence of small anatomical geometric complexities such as the orbital cavities to positively influence resolution using elastic waves. We believe that it will be useful, in future studies, to explore the performance of our algorithms in other frequency ranges and for other biological models (e.g., echolocating species such as dolphins or bats).

As noted by Parseihian et al. [START_REF] Parseihian | Reaching nearby sources: comparison between real and virtual sound and visual targets[END_REF], very few studies in psychoacoustics have explored human sound lo-It appears to us that further experimental work is 396 needed to more robustly evaluate how well humans 397 localize nearby sources and if our findings can be re-398 lated to psychoacoustic studies in the near-field.

Figure 1 :

 1 Figure1: Sketch of the experimental setup in the horizontal plane: A loudspeaker is connected to a source generator (PC) and emits a chirp signal at each angle ϕ ranging from -50 • to 50 • along a half circle at various distances to the skull. The resulting vibration of the skull is recorded through two passive sensors glued to the hypothetical ear locations. They are connected to the signal acquisition system, consisting of a sound card connected to a PC.

164(Figure 2

 2 Figure 2: a) Exemplary waveform of a recorded signal at one of the sensors. b) Frequency spectrum of the same signal.

1 -Figure 3 :Figure 4

 134 Figure3: Normalized focusing functions along the curvilinear abscissa for sources in front of the center of the skull (ϕ = 0 • ) and at different distances to the skull. The distance of the measurement points to the skull decreases from 40 cm, down to 20 cm, 12 cm and 5 cm (different curves). There is a clear trend of increasing resolution (decreasing -3 dB width of the curves) with decreasing distance.

Figure 4 :

 4 Figure 4: -3 dB width values of the focusing functions for sources at different distances to the skull (x-axis) and maximum frequencies f max of the signal. The slope of each linear fit, which corresponds to the apparent aperture A in Equation 4, is approximately 10 cm for all curves.

Figure 5 :

 5 Figure 5: Angular variations of resolution for different source distances. Top: Variation in apparent aperture for different source distances. Maxima are at -20 • and 15• whereas the values decrease for source positions close to the center and further away from the center. Bottom: Variation in -3 db widths for different source distances. Super-resolution is accomplished throughout all angles at a distance of 5 cm and for certain angles at a distance of 12 cm. Highest resolution (smallest -3 dB width) is accomplished for source positions directly in front of the orbital cavities. This effect is (relatively) enhanced the closer the source to the skull.

327- 3

 3 dB widths are smaller than half a wavelength 328 (super-resolution) throughout all azimuths at a dis-329 tance between source and skull of 5 cm, down to 330 λ min /15 (i.e. for ϕ =-20 • and 15 • ). This shows that 331 the skull-shaped antenna enables sub-wavelength fo-332 cusing of near-field sources and, furthermore, anatom-333 ical details of the skull may give rise to differences in 334 resolution at certain positions due to the presence of 335 evanescent waves. They can be described as a non-336 propagative spatial fluctuation field that decreases ex-337


Magdalena Oryaëlle Chevrel
email: oryaelle.chevrel@gmail.com

Andrew J L Harris

Mike R James

Laura Calabrò

Lucia Gurioli

Harry Pinkerton

The viscosity of pāhoehoe lava: in situ syn-eruptive measurements from Kilauea, Hawaii

Keywords: lava, rheology, rotational viscometer, bubbles, crystals, basalt

Viscosity is one of the most important physical properties controlling lava flow dynamics. Usually, viscosity is measured in the laboratory where key parameters can be controlled but can never reproduce the natural environment and original state of the lava in terms of crystal and bubble contents, dissolved volatiles, and oxygen fugacity. The most promising approach for quantifying the rheology of molten lava in its natural state is therefore to carry out direct field measurements by inserting a viscometer into the lava while it is flowing. Such in-situ syn-eruptive viscosity measurements are notoriously difficult to perform due to the lack of appropriate instrumentation and the difficulty of working on or near an active lava flow. In the field, rotational viscometer measurements are of particular value as they have the potential to measure the properties of the flow interior rather than an integration of the viscosity of the viscoelastic crust + flow interior. To our knowledge only one field rotational viscometer is available, but logistical constraints have meant that it has not been used for 20 years. Here, we describe new viscosity measurements made using the refurbished version of this custom-built rotational viscometer, as performed on active pāhoehoe lobes from the 61G lava flow of Kilauea's Pu'u 'Ō'ō eruption in 2016. We successfully measured a viscosity of ~380 Pa s at strain-rates between 1.6 and 5 s -1 and at 1144 °C. Additionally, synchronous lava sampling allowed us to provide detailed textural and chemical characterization of quenched samples. Application of current physico-chemical models based on this characterization (16±4 vol.% crystals; 50±6 vol.% vesicles), gave viscosity estimates that were approximately compatible with the measured values, highlighting the sensitivity of model-based viscosity estimates on the effect of deformable bubbles. Our measurements also agree on the range of viscosities in comparison to previous field experiments on Hawaiian lavas. Conversely, direct comparison with sub-liquidus rheological laboratory measurements on natural lavas was unsuccessful because recreating field conditions (in particular volatile and bubble content) is so far inaccessible in the laboratory. Our work shows the value of field rotational viscometry fullyintegrated with sample characterization to quantify three-phase lava viscosity. Finally, this work suggests the need for the development of a more versatile instrument capable of recording precise measurements at low torque and low strain rate, and with synchronous temperature measurements.

Introduction

Understanding and describing the complexity of lava flow behavior is a major challenge and a long-term objective in modern volcanology. Modeling of lava flows leads to significant improvements in hazard assessment, as well as contributions to our understanding of volcanic activity and history on other planets. Rheology, which is directly linked to the intrinsic chemical and physical properties of the magma (chemical composition, oxygen fugacity, volatile content, temperature, and shape and size of bubbles and crystals) is a key influence on the transport of magma from its source to the surface (e.g., [START_REF] Dingwell | Volcanic Dilemma Flow or Blow[END_REF]. As lava flows, cooling and degassing trigger crystal and bubble growth and interstitial melt differentiation, driving constant evolution of the chemical and physical properties (e.g., [START_REF] Crisp | Crystallization History of the 1984 Mauna Loa Lava Flow[END_REF]Applegarth et al. 2013). This continuous transformation leads to variations in rheological properties (viscosity, yield strength, strain-rate dependency) which directly impact the dynamics of lava flow emplacement (e.g., Pinkerton and Norton 1995;[START_REF] Harris | FLOWGO: A Kinematic Thermo-Rheological Model for Lava Flowing in a Channel[END_REF].

Several lava flow emplacement models exist in the literature (see [START_REF] Cordonnier | Benchmarking Lava-Flow Models[END_REF] for comparison between the different models). Most of these models require parameterization of the lava physical properties (such as viscosity) and of the eruption conditions (such as eruption temperature or effusion rate), all of which are challenging to measure. Lava rheology can be determined using various methods. Following the pioneering study of [START_REF] Nichols | Viscosity of Lava[END_REF], lava rheological parameters can be estimated from flow velocity, underlying slope and channel dimensions assuming Newtonian or Bingham rheology. The advantage of this method is that channel velocity is a parameter that can be measured at any type of flow and provides the flow-scale 'bulk rheology' (e.g., [START_REF] Moore | Preliminary Estimates of the Rheological Properties of 1984 Mauna Loa Lava[END_REF][START_REF] James | Image-Based Measurement of Flux Variation in Distal Regions of Active Lava Flows[END_REF]). However, this method cannot quantify the heterogeneity of rheological properties within and across a flow.

Another approach is to perform laboratory measurements using re-melted natural samples at super-and sub-liquidus conditions. Recent examples include [START_REF] Ishibashi | Non-Newtonian Behavior of Plagioclase-Bearing Basaltic Magma: Subliquidus Viscosity Measurement of the 1707 Basalt of Fuji Volcano, Japan[END_REF], [START_REF] Vona | The Rheology of Crystal-Bearing Basaltic Magmas from Stromboli and Etna[END_REF], [START_REF] Vetere | Viscosity Changes during Crystallization of a Shoshonitic Magma: New Insights on Lava Flow Emplacement[END_REF], [START_REF] Chevrel | Viscosity Measurements of Crystallizing Andesite from Tungurahua Volcano (Ecuador)[END_REF], [START_REF] Sehlke | Pahoehoe to `a`a Transition of Hawaiian Lavas: An Experimental Study[END_REF], and [START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF]. These measurements are usually performed using a concentric cylinder apparatus, and the full rheological curve may be constructed for a given crystallinity and temperature. Recent experiments have explored the effect of cooling rate on crystallinity and the rheological response of the lava [START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF][START_REF] Kolzenburg | The Rheological Evolution of the 2014/2015 Eruption at Holuhraun, Central Iceland[END_REF]. Laboratory measurements have the advantage that they are made in an environment that allows key parameters to be controlled and precise rheological measurements to be made. However, they are unable to replicate the range of crystal and bubble contents, dissolved volatiles, and oxygen fugacity of natural lavas.

When laboratory experiments cannot be carried out, or when an active flow has not been observed, a common practice is to estimate viscosity from a "petrologic" approach. This involves applying empirical physico-chemical models based on previous laboratory experiments to estimate viscosity from groundmass chemical composition and crystal content (e.g., [START_REF] Pinkerton | Methods of Determining the Rheological Properties of Magmas at Sub-Liquidus Temperatures[END_REF]. This approach is clearly the best long-term method of providing rheological data for modeling volcanic processes because, being based on sample texture characterization, it can be applicable to present, past or extraterrestrial deposits. The petrologic approach may be applied to quenched samples from an active flow, where the groundmass crystallization will evolve as a function of cooling through time and distance from the vent (e.g., [START_REF] Crisp | Crystallization History of the 1984 Mauna Loa Lava Flow[END_REF][START_REF] Riker | The Length of Channelised Lava Flows: Insight from the 1859 Eruption of Mauna Loa Volcano, Hawaii[END_REF]. The effect of bubbles (term used when describing the molten lava) may also be taken into account from estimation of the vesicle (term used when describing the rock) content, which may also vary down flow [START_REF] Robert | Textural and Rheological Evolution of Basalt Flowing down a Lava Channel[END_REF]. Another use of this method is to determine the viscosity from the crystallization sequence as calculated at thermodynamic equilibrium (e.g., [START_REF] Harris | FLOWGO: A Kinematic Thermo-Rheological Model for Lava Flowing in a Channel[END_REF][START_REF] Riker | The Length of Channelised Lava Flows: Insight from the 1859 Eruption of Mauna Loa Volcano, Hawaii[END_REF][START_REF] Chevrel | Lava Flow Rheology: A Comparison of Morphological and Petrological Methods[END_REF]. However, care must be taken when assuming thermodynamic equilibrium because lava flows are undercooled systems whose properties (crystallization and hence rheology) are cooling-rate dependent [START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF]. During flow, lava does not reach the crystal/melt equilibrium and this disequilibrium leads to a delay in crystallization that must be accounted for when calculating lava viscosity down flow.

Currently, the only way to determine the complex rheological behavior of lava in its natural state is to measure it directly in the field by inserting a viscometer into the molten lava while it is flowing. Such in-situ viscosity measurements are notoriously challenging to perform due to the difficulty of approaching an active lava flow, and the lack of appropriate instrumentation. A small number of viscosity measurements have been made using either penetrometers [START_REF] Einarsson | The Flowing Lava. Studies of Its Main Physical and Chemical Properties[END_REF][START_REF] Gauthier | Field and Laboratory Studies of the Rheology of Mount Etna Lava[END_REF][START_REF] Panov | Mechanical Properties of Lavas Extruded in the 1983 Predskazannyi Eruption (Klyuchevskoy Volcano)[END_REF][START_REF] Pinkerton | Field Measurements of the Rheology of Lava[END_REF]Belousov and Belousova 2018) or rotational viscometers (Pinkerton et al. 1995;[START_REF] Pinkerton | Rheological and Related Properties of Lavas[END_REF]). Among these studies, measurements performed using simple lava penetrometers (a metal rod pushed into the lava) produce results that are influenced by the outer viscoelastic layer. Measurements start as soon as the sensor is inserted into the lava so that the viscosity of the viscoelastic crust + flow interior is integrated. However, all rotational viscometers, and the custom-build penetrometer used by [START_REF] Pinkerton | Field Measurements of the Rheology of Lava[END_REF], have the potential to measure the properties of the flow interior because measurements start once the isothermal interior is reached. The viscous core of large 'a'ā flows remains inaccessible for field viscometry, due to the problems of penetrating the overlying breccia layer. Low viscosity lava flows, such as pāhoehoe type and small channelized flows are the only targets so far accessible by field viscometry.

To our knowledge, viscometers developed by the University of Lancaster in the 1980's and 1990's [START_REF] Pinkerton | Rheological and Related Properties of Lavas[END_REF]Pinkerton et al. 1995a,b) are the only portable rotational viscometers to have been used in situ on active lava flows, with the last measurements performed more than 20 years ago [START_REF] Pinkerton | Rheological and Related Properties of Lavas[END_REF]Pinkerton et al. 1995b). Since then no further work has been presented due to lack of suitable conditions and personnel (and funds) dedicated to the project. Conversely, over the past decades, extensive laboratory measurements on natural or analogue mixtures have greatly improved our understanding of the chemical control of the melt phase on viscosity [START_REF] Giordano | Viscosity of Magmatic Liquids: A Model[END_REF], as well as the rheology of two-phase mixture [START_REF] Mader | The Rheology of Two-Phase Magmas: A Review and Analysis[END_REF]. As a result, although many physicochemical based models now exist to estimate the rheology of magmas and lavas, no field measurements have been performed to test whether, and how, they can be applied to natural lava flows. Additional field validation is therefore required to provide confidence in the petrologic approach.

In this work, we restored the unique field rotational viscometer of Pinkerton et al. (1995b) and performed viscosity measurements on pāhoehoe lobes at Kilauea, Hawaii. The aims of this paper are to present this refurbished instrument and our new field measurements.

Our viscometry results are described in detail, and complemented by textural and chemical analyses performed on quenched samples, which we use to implement, and compare with, a petrologic approach.

Methods

Viscometry

Description of the field viscometer

We restored the rotational viscometer used by Pinkerton et al. (1995b) for field measurements (Fig. 1). This instrument consists of a 24-V DC variable speed Bosch drill that drives a vane, which is inserted into the lava. The vane geometry is unlike many laboratory instruments, where the immersed spindle is often a cylinder. Instead it was designed for field use to minimize disturbance of the fluid during immersion whilst minimizing slippage between the spindle and lava. The vane is attached to the rotating inner shaft which is protected with a fixed outer tube equipped by bearing assemblies (containing graphite rings to minimize friction). This helps to maintain alignment and low-friction rotation of the inner shaft (Fig. 1a), which is driven by the drill via a 15:1 reduction gearbox, a torque limiter (2 Nm), and a torque and rotation rate sensor (TORQSENSE E300 Rayleigh Wave Transducer from Sensor Technology Ltd). The sensor is linked to a stand-alone E302 interface/readout connected to a laptop via USB, and is a different version to the one in Pinkerton et al. (1995b). In 2016, the torque sensor and transducer were calibrated by Sensor Technology and small adjustments in communication were made so that the computer simultaneously recorded the rotation rate and torque every 0.05 s through the software provided (TORQVIEW2). At rest, we recorded a sensor sensitivity of ±0.02 Nm. When a rotation was applied with the instrument mounted with the vane in air, an oscillation of ±0.05 Nm was recorded, although efforts were made to reduce this by improving the alignment of the rotating rod between the motor and the vane. New shear vanes were made using the same type of high temperature-resistant stainless steel (BS 321 S20 or Z6CNT18.10) as used for the inner shaft and outer tube. Two vane sizes were employed, both composed of four paddles (2 mm thick) welded at 90° onto a 1 cm diameter rod, so that the vane diameter across the fins is 40 mm for the 60 mm-long vane (hereafter S60/20), and 50 mm for the 80 mm-long vane (hereafter S80/25). The total length of the assembled apparatus is 2.7 m, has a weight of about 15 kg, and can be dismantled into four pieces for transport.

Data processing to extract viscosity

The theory employed is that of wide-gap concentric cylinder viscometry. Here, the torque is converted into shear stress and rotational speed into strain-rate using the spindle geometry via the Couette theory for Newtonian or non-Newtonian fluids [START_REF] Spera | Rheology of Melts and Magmatic Suspensions I. Design and Calibration of a Concentric Cylinder Viscometer with Application to Rhyolitic Magma[END_REF]Pinkerton and Norton 1995;[START_REF] Stein | New High-Temperature Rotational Rheometer for Silicate Melts, Magmatic Suspensions, and Emulsions[END_REF][START_REF] Vona | The Rheology of Crystal-Bearing Basaltic Magmas from Stromboli and Etna[END_REF]. Material trapped between the paddles will move as the spindle rotates and therefore a virtual cylinder of sample material is used for the calculation. The shear stress is calculated via:

! = ! 2!ℎ! ! !
where M is the torque recorded by the torque sensor, h is vane length and R i is the equivalent radius of the rotating vane. Strain rate is obtained from the applied angular velocity via [START_REF] Stein | New High-Temperature Rotational Rheometer for Silicate Melts, Magmatic Suspensions, and Emulsions[END_REF]:

! = 2Ω ! 1 - ! ! ! ! ! ! (2)
where Ω is angular velocity (in rad/s), n is the flow index, and R o is the radius of the outer cylinder. In the case of field measurements, the vane is effectively immersed into an unconstrained medium (i.e., R o approaches infinity). Now, following Barnes (1989;see supplement), the calculation of strain rate (for a range of 0.1 to 10 s -1 ) can be reduced to:

When several angular speeds are applied, n is usually obtained by calculating the slope of the measured ln(τ) vs. ln(Ω) relation. The flow curves (i.e., the graph of strain rate versus shear stress) can then be established to determine the fluid's rheological model. This is Newtonian when viscosity is proportional to strain rate (n = 1):

or non-Newtonian when viscosity varies with strain rate (n ≠ 1) as, for example, in a power law model:

Here K is the flow consistency, which corresponds to viscosity at unit strain rate.

In practice, when measuring rheology, if one or only a few angular velocities are applied, the full flow curve cannot be reliably established so that n cannot be accurately determined. A common practice is therefore to calculate the apparent viscosity that is equivalent to the Newtonian viscosity. This can be determined at each applied strain rate by assuming n = 1.

Testing the field viscometer in the laboratory

The instrument was tested in the laboratory using the viscosity reference standard N190000 from Cannon Instrument Company®. This has been certified to be 33.06, 143.8, 523.3 and 843 Pa s at 60, 40, 25 and 20 °C, respectively. These certified values were used to define the exponential relationship: η = 3990.7e -0.081T (R 2 =0.99), allowing us to estimate the viscosity of the N190000 standard at experimental temperatures. To perform the measurements, 1.5 liters of N190000 viscosity standard were poured into a plastic container, which was 12 cm in diameter and 15 cm deep. Before each experiment, the fluid was placed overnight in an oven at the desired temperature to allow equilibration and removal of all air bubbles. Prior to each calibration, the container was removed from the oven and the temperature of the fluid measured using a K-type thermocouple located in the liquid beside the vane. Temperatures across the container varied by less than 2 °C and were constant near the vane during the experiment, with a maximum variation of ±0.6 °C. The instrument was oriented at an angle of 45° and the container was inclined so that the vane was fully immersed into the liquid. Both vanes were tested at several temperatures between 20 and 40 °C and over rotation speeds of 2

! = 2Ω ! (3) ! = ! ! (4) ! = !! ! (5)
to 20 rpm, resulting in strain rates adjacent to the vane of 0.6 to 5 s -1 (calculated using Eq. 2). Higher strain rates were not applied because of shear thinning effects that will bias the measurements (Cannon Instrument Company, personal communication, October 2016). The average torque was obtained over 8 to 60 s of stable reading at different strain rates at each temperature (Fig. 2 and supplement). The ln(τ) vs. ln(Ω) trends (Fig. 2b) had slopes of n = 1±0.09 as expected for Newtonian behavior of the N190000 viscosity standard. Apparent deviation from Newtonian behavior was observed when torque was less than 0.1 Nm (equivalent to stresses of 318 and 663 Pa for vanes S80/25 and S60/20, respectively). This is attributed to limits in the torque transducer capability at low torque values, and therefore only higher torque measurements were considered. Viscosity was then calculated following wide gap theory (Eqs. 1, 2 and 4) for Newtonian behavior. Less than 5 % of error and less than 50 Pa s of deviation from the standard values was found (Fig. 2c).

Operating the field viscometer in the field

To perform a measurement in the field, two people held the instrument. The operator closest to the lava inserted the vane into the active lobe and held it in position in the middle of the fluid core taking care not to ground the vane. The second person controlled the drill speed and monitored the data being collected on the computer. When the targeted lava lobe was surrounded by other hot surfaces, the computer would be monitored by a third person located a couple of meters away in cold air and who would direct operations. At the same time, temperature measurements were made with a K-type thermocouple (Fig. 1b). This was done only for the first measurement because, after this, the third person was needed to perform other measurements. Before insertion, the vane was pre-heated by being placed close to the lobe surface and exposed to the radiant heat. However, this pre-heating procedure did not last more than few seconds and could not always be performed because the active lobe cooled too rapidly so that an impenetrable outer crust formed, or because operators could not handle the radiant heat themselves. The instrument was then rapidly inserted into the hottest and most fluid part of the active lobe and the measurement started. As the lava advanced and its surface began to solidify, the instrument was slowly pulled backwards by the operator. Each measurement lasted less than one minute, mainly due to lava cooling around the rod to prevent rotation. Radiant heat was also a limiting factor for the operators, despite the use of heat protective suits and hoods (Fig. 1). Although we aimed to reach a stable reading over 20 s (as in the laboratory), it was not always possible because the surface cooled too quickly and imprisoned the spindle. To record a range of strain rates, the angular velocity was manually increased ("up path") and then decreased ("down path") by changing the drill speed. However, practicalities meant that it was not always possible to cover the full up and down path ranges. To estimate the apparent viscosity, the torque recorded during constant drill speed over at least 3 to 23 s was then averaged and the standard deviation from the mean used for error bars (see example in Fig. 3a and raw data in supplement).

Lava sampling and textural analyses

Samples were collected after each viscosity measurement by quenching the lava attached to the vane in a water bucket. Another set of samples was collected using a stainless steel tube inserted into an active pāhoehoe lobe. The tube was inserted into the lobe front and withdrawn from the lava before quenching with water, effectively coring the tube. Additionally, one control sample was collected from a pāhoehoe lobe emplaced a few weeks earlier. Bulk rock major element analysis of one sample was carried out with Inductively Coupled Plasma -Atomic Emission Spectroscopy at the Laboratoire Magmas et Volcans (LMV, Université Clermont-Auvergne, France). Thin sections were made from all samples, and textural and chemical analyses were carried out at LMV. Olivine and glass chemical compositions were obtained from around 10-20 analyses for each thin section via the electron microprobe (CAMECA SX 100), operating at 15 kV and with a focused beam of 15 nA for olivine and a 20 µm defocused beam of 8 nA for glass. Total water content was determined from double polished samples by Fourier Transform Infrared Spectroscopy, using a Bruker Vertex 70 spectrometer coupled with a Hyperion microscope system, following the methods given in Mercier et al. (2010, see supplement). Crystallinity, vesicularity and vesicle size distribution were estimated from analyses of the thin sections through FOAMS following the methods described in [START_REF] Shea | Textural Studies of Vesicles in Volcanic Rocks: An Integrated Methodology[END_REF]. For each thin section, high magnification images were acquired with scanning electron microscopy (SEM) in BSE (backscattered electron imagining mode) at ×25 magnification (six images in total), and ×100 (fifteen images in total). The smallest object measured was 10 pixels across (equivalent to 0.012 mm).

A mean value (averaged from two to four measurements) for the vesicle-free rock density (DRE density) was determined by powdering the rock and measuring the volumes of known masses using an Accupyc 1340 Helium Pycnometer. The bulk density of the samples was measured using a pycnometer (Micromeritics Geopyc 1360 envelope density analyzer). This instrument calculates the envelope density of the sample, which is the mass divided by its encompassing volume, where the encompassing volume includes pores and small cavities. The measurement technique consists of measuring the difference in the volume of a quasifluid medium (DryFlo TM , composed of tiny, rigid spheres) with and without the sample embedded in the medium. Prior to immersion samples were wrapped with parafilm® to preserve external irregularities but to avoid the medium entering the porous sample. The density-derived vesicularity was then calculated using the relationship given by [START_REF] Houghton | A Vesicularity Index for Pyroclastic Deposits[END_REF].

Results

The targeted lava flow was the 61 G flow of Kilauea's Pu'u 'Ō'ō eruption (Hawaii) which had been active since May 2016. At the time of the field campaign, the lava was mainly flowing in tubes from the vent down to the ocean entry, but on November 23, 2016, we were able to deploy the instrument on several pāhoehoe lobes at N19.349537, W155.048793. These had broken out from the tube 7.2 km from the vent overnight. The lava flow field consisted exclusively of S-type pāhoehoe lava with a bubble-rich silvery surface. Each observed lobe was 10 to 30 cm thick upon emplacement and was active for about three to five minutes before stalling, inflating and erupting another lobe from its base. Core temperature was measured at 1144 °C with a K-Type thermocouple inserted to 5 to 8 cm into the lobe (Fig. 1b).

Viscometry

Viscosity measurements were successfully performed on seven distinct lobes. For each run, between one and four strain rates associated with stable torque readings could be collected before the spindle stopped rotating due to lobe cooling (the example of run 12 is given in Fig. 3a) or due to operator discomfort because of having to stand on an active surface or strong radiative heat. For four runs (R2, R9, R10, R11), data could also be collected during the down path but showed a marked hysteresis when compared with the up path, resulting in a higher apparent viscosity (see supplement). A similar viscosity increase was also recorded for some measurements within the up path phase (for example R2, R8, R9), which could be interpreted as representing shear thickening. However, shear thickening has not been observed in previous lava rheology measurements, so the effect is more likely to be due to rapid cooling of the lobe, or to the accumulation of cooled lava around the vane. Both would have significantly increased the effective diameter of the spindle. Considering data that were not affected by this phenomenon, rotation velocities ranged from 7.6 to 23.8 rpm, corresponding to strain rates between 1.6 and 5 s -1 at the vane. Torque was recorded between 0.11 and 0.55 N m, reflecting shear stresses of between 755 and 2132 Pa (Table 1). The apparent viscosity at each strain rate was determined to be between 261 and 483 Pa s. Unfortunately, the limited number of measurements possible at each lobe did not allow us to construct the full rheological flow curve with confidence. We also note that no yield strength was measured. Gathering all runs into a single plot (Fig. 3b) shows that the lava had a Newtonian viscosity (Eq. 4) of 379 Pa s (R 2 =0.83). Treating the data as a power-law rheology (Eq. 5) gave a flow index of 0.88 and a consistency of 424 Pa s (R 2 = 0.79).

Sample chemical analyses and textural characterization

Bulk rock analysis revealed 50.8 wt. % SiO 2 and a Mg# of 0.38 (Table 2), with a DRE density of 2847 kg m -3 . The interstitial glass chemical composition was identical and homogenous in all samples and slightly more evolved than the bulk rock, with a water content of 0.077±0.015 wt.% (Table 2). All samples had a porphyritic texture containing phenocrysts mainly of olivine and rarely of plagioclase, along with microlites (<0.1 mm) of olivine and plagioclase, as well as minor pyroxenes and minor oxides (usually as inclusions in olivine phenocrysts). These were embedded within a glassy matrix (Fig. 4). Olivine crystals had euhedral and subeuhedral shapes and in some samples (CN, C1a and R9), they were organized in glomerocrysts with plagioclase (Fig. 4c). Plagioclase occurs mostly as microlites (in some samples all plagioclase were as microlites) and pyroxene was the least abundant phase and occured only as microlites. The aspect ratio (length vs. width) of the mafic minerals was 1.8 and 3.1 for plagioclase. The vesicle-free crystal fraction (phenocrysts + microlites) varied from 0.12 to 0.21 (Table 3). Overall, the similarity in crystal shape and content between the control sample and the quenched samples indicated that our sampling technique did not affect the texture of the lava. Therefore, the quenched samples could be used to estimate the effect of crystals on viscosity.

In terms of vesicularity, the quenched samples showed some differences with the control sample. The control sample had a bulk density of 1441 kg m -3 , which lead to a density-derived vesicularity of 0.49. Image analysis of the control sample revealed a 2D vesicularity of 0.52 (Table 3). In contrast, quenched samples had lower densities (averaged at 1161±171 kg m -3 ) and therefore higher density-derived vesicularities (averaged at 0.60±0.06; Table 3). The difference in vesicularity between the control sample and the quenched samples was due to the presence of isolated vesicles larger than 8 mm in the latter (see R12 in Fig. 4). Although not always observed in thin section, these large vesicles were always present in the quenched samples (as observed in hand specimen and as inferred from the higher density-derived vesicularity). The formation of these large vesicles was likely due to the ingestion of air in the molten lava during sampling, either when the tube was inserted or when the vane was withdrawn. To evaluate the bubble fraction in the flowing lava and during viscosity measurements, the vesicularity of the quenched samples therefore needed to be corrected for these large vesicles that represent artifacts introduced by the sampling process. Vesicle size distributions revealed that, for all the samples, there was a common population of spherical (when single) to convoluted (when two or more bubbles had coalesced) vesicles that range between 0.1 and 5 mm in size (Fig. 5). The average vesicle fraction obtained from the images (Table 3), which is that used here to estimate the effect of bubbles on viscosity, was 0.50±0.06. All samples also contained a population of irregular micro-vesicles smaller than 0.03 mm. These were probably diktytaxitic voids related to the crystallization of microlites. We note that samples quenched after the viscosity measurement (R9, R10 and R12) also had a vesicle population between 0.03 to 0.1 mm that was not observed in either the control sample or the quenched sample from the sampling tube (Fig. 5). This population may have formed because of shearing during viscometry that led to disruption of larger bubbles to form smaller bubbles [START_REF] Stein | Rheology and Microstructure of Magmatic Emulsions: Theory and Experiments[END_REF]).

Comparison with viscosity estimated from other techniques

Petrologic approach using textural and chemical analyses

The parameters used for viscosity estimations via the petrologic approach are given in Table 4. Using the glass composition, and including the H 2 O content (Table 2) and the measured temperature (1144 °C), the viscosity of the melt phase was calculated using the model of [START_REF] Giordano | Viscosity of Magmatic Liquids: A Model[END_REF] to be 330 Pa s. Note that a variation of ±1 °C affected the melt viscosity by 2 %. Considering the average vesicle-free crystal fraction, the mixture (melt + crystals) viscosity was estimated via the method described by [START_REF] Mader | The Rheology of Two-Phase Magmas: A Review and Analysis[END_REF] for rough particles with an aspect ratio of 2.4, this being the average value calculated for all the crystals. The resulting bubble-free mixture viscosity is 699 Pa s. Note that a variation of ±1 vol. % in crystalinity changed the viscosity by 6 %. The effect of bubbles on the mixture viscosity depends on their ability to deform. We calculated the bubble capillary number (Ca) from the applied shear rates during the viscosity measurements (3.5 s -1 on average), with a modal bubble radius of 0.34 mm. Using the melt viscosity as calculated above and a bubble-liquid interfacial tension of 0.3 Nm -1 [START_REF] Murase | Properties of Some Common Igneous Rocks and Their Melts at High Temperatures[END_REF], we obtained a Ca of 1.3. Following [START_REF] Llewellin | Bubble Suspension Rheology and Implications for Conduit Flow[END_REF], 50 vol. % of deformable bubbles would lower the viscosity by a factor of 0.31, to give a bulk viscosity of the three-phase mixture (melt+crystals+bubbles) of 220 Pa s. Accounting for the variability in crystal (±4 vol. %), bubble (±6 vol. %) and water (±0.015 wt. %) content this value varied by ±130 Pa s. Thus, when accounting for deformable bubbles, the petrologic model provided a reasonable approximation to the field-measured viscosity. In contrast, neglecting the effect of bubbles, or considering bubbles as non-deformable solid particles, would lead to a significant overestimation in viscosity (this being >10 3 Pa s). We note that the differences in viscometermeasured viscosity between each lobe did not exactly correlate with the model-based viscosity estimation made from the crystal and bubble content of the quenched samples (Fig. 3). This can either be explained by our field viscosity measurements needing better accuracy, or because the rheological models used in the petrologic approach were inappropriate.

Laboratory results from previous studies

Laboratory viscometry of crystallizing Hawaiian lavas at subliquidus temperatures were performed by [START_REF] Shaw | Rheology of Basalt in the Melting Range[END_REF], [START_REF] Ryerson | Rheology of Subliquidus Magmas: Picritic Compositions[END_REF], and [START_REF] Sehlke | Pahoehoe to `a`a Transition of Hawaiian Lavas: An Experimental Study[END_REF]. The main difference between these experiments and field measurements is the absence of bubbles in the laboratory experiments. Another important difference is that crystallization in the laboratory always takes place at a higher temperature than in nature and, consequently, viscosity will always be greater at any given temperature in the laboratory experiments. This is mainly due to three reasons. First, laboratory experiments are usually performed under atmospheric conditions at 1 atm with an oxidized oxygen fugacity and a completely degassed melt. Second, constant stirring of the melt during viscosity measurement enhances crystallization [START_REF] Kouchi | Effect of Stirring on Crystallization Kinetics of Basalt: Texture and Element Partitioning[END_REF][START_REF] Vona | The Effects of Undercooling and Deformation Rates on the Crystallization Kinetics of Stromboli and Etna Basalts[END_REF][START_REF] Chevrel | Viscosity Measurements of Crystallizing Andesite from Tungurahua Volcano (Ecuador)[END_REF]. Third, most laboratory experiments apply isothermal conditions and are at thermodynamic equilibrium. This results in higher crystal contents than under steady cooling (disequilibrium) as in nature [START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF]. For example, at 1144 °C (the lava temperature measured in the field), experiments by [START_REF] Sehlke | Pahoehoe to `a`a Transition of Hawaiian Lavas: An Experimental Study[END_REF] would predict a crystal content of almost 100% and a viscosity of more than 10 5 Pa s. Comparing the experiments of [START_REF] Shaw | Rheology of Basalt in the Melting Range[END_REF] and [START_REF] Ryerson | Rheology of Subliquidus Magmas: Picritic Compositions[END_REF] at this same temperature, viscosity would be predicted to be almost the same as we measured but with a much higher crystal content [25 and > 45 vol. % for [START_REF] Shaw | Rheology of Basalt in the Melting Range[END_REF] and [START_REF] Ryerson | Rheology of Subliquidus Magmas: Picritic Compositions[END_REF], respectively].

Following [START_REF] Ryerson | Rheology of Subliquidus Magmas: Picritic Compositions[END_REF] the lava should have a yield strength of 28 Pa and 125 Pa for 15 and 25 vol. % crystals, respectively, while following [START_REF] Sehlke | Pahoehoe to `a`a Transition of Hawaiian Lavas: An Experimental Study[END_REF], the yield strength should be lower. However, whether a crystal content as low as 15 to 25 vol.% would build a crystal framework and develop a yield strength is a matter of ongoing debate [START_REF] Saar | Numerical Models of the Onset of Yield Strength in Crystal-Melt Suspensions[END_REF]. Any yield strength present during our experiments was smaller than the level of detectability of the instrument which means that -if present -the yield strength is less than 300 Pa. Lower strain rates and a more sensitive torque sensor are required to improve measurement sensitivity. Additionally, better constraints on estimating the yield strength from the petrologic approach and that considers crystals and bubbles are needed (see discussion in [START_REF] Chevrel | Lava Flow Rheology: A Comparison of Morphological and Petrological Methods[END_REF]Castruccio et al. 2014).

Previous field-based viscometry measurements

Two previous in situ field viscosity measurements were undertaken on Hawaiian lavas. The first measurements were made at Makaopuhi lava lake by [START_REF] Shaw | The Viscosity of Basaltic Magma: An Analysis of Field Measurements in Makaopuhi Lava Lake, Hawaii[END_REF] using a vertical rotational viscometer. They recorded a viscosity of 650-720 Pa s for a unit strain rate and yield strengths of 70 and 120 Pa. The collected samples and in situ temperature measurements revealed <5 % vesicles and a crystallinity of 25% at 1130±5 °C. Although our samples had a higher silica content, our viscosity measurements were slightly lower. This is explained by a combination of lower crystal fraction, higher content of deformable bubbles and higher temperature. The other set of measurements was performed by Pinkerton et al. (1995) on three 0.2 to 0.5 m thick pāhoehoe lobes erupted in September 1994 at Kilauea using the same rotational viscometer as used in this study but equipped with a more sensitive torque sensor. Maximum measured temperatures were 1146 °C and apparent viscosities at unit strain rate were in the range 234-548 Pa s. The lowest viscosity lobe was also the closest to Newtonian behavior and was best characterized by a power law model of the form ! = 234! !.!! . The highest viscosity lobe departed significantly from Newtonian behavior and fitted a power law model with the form ! = 548! !.!" . The non-Newtonian behavior measured in 1994 might reflect different vesicularities or crystallinities compared with those of the 2016 lavas, however samples collected in 1994 were not analyzed. Alternatively, it could be the result of operating the instrument at lower rotational speeds and hence strain rates (< 1 s -1 ) in 1994. This permitted construction of a complete rheological curve (Fig. 3). Although some differences were observed between the 1994 and 2016 measurements, the results of both studies agreed on the range in apparent viscosities at unit strain rate (Fig. 3).

Discussion

We have compared our field viscosity measurements with estimations of lava viscosity from the petrologic approach and from previous laboratory and field measurements. However, the measured values represent only a "snap shot" of the lava interior behavior over a short time scale and at a single temperature. Indeed, pāhoehoe lobes cool very quickly (e.g., [START_REF] Hon | Emplacement and Inflation of Pahoehoe Sheet Flows: Observations and Measurements of Active Lava Flows on Kilauea Volcano, Hawaii[END_REF]Ball et al. 2008;[START_REF] Gottsmann | Thermal History of Hawaiian Pāhoehoe Lava Crusts at the Glass Transition: Implications for Flow Rheology and Emplacement[END_REF]). Ball et al. (2008) showed that small stationary lobes (without advection of hot lava) and non-stationary lobes cool to 800 °C in 25 s and 50 s, respectively. Thus, lobe cooling will lead to a rapid increase in viscosity that will quite quickly impede viscosity measurement. This was observed during our measurements when the vane became stuck or when viscosity measurements dramatically increased, sometimes after only 10 s. In general, we observed that viscosity always increased with time; which likely reflected the cooling effect of the whole lobe. The increase of viscosity with time, could also be due to accumulation of lava around the vane due to the chilling effect of the vane being inserted, which increased the apparent vane diameter. To discriminate between lobe cooling or accumulation effects, further field viscosity measurements need to be performed while continuously recording the temperature at the vane. For this, the rotating vane should be equipped with a thermocouple, as has been recently achieved in the laboratory [START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF]. This type of measurement will allow tracking of the temperature, and therefore the cooling rate -viscosity relationship, which is a key to understanding the thermo-rheological control of lava flow emplacement [START_REF] Giordano | Thermo-Rheological Magma Control on the Impact of Highly Fluid Lava Flows at Mt. Nyiragongo[END_REF][START_REF] Kolzenburg | The Rheological Evolution of the 2014/2015 Eruption at Holuhraun, Central Iceland[END_REF][START_REF] Kolzenburg | In Situ Thermal Characterization of Cooling/crystallizing Lavas during Rheology Measurements and Implications for Lava Flow Emplacement[END_REF]. Another problem that was faced during our measurements was the cooling of the lobe crust around the rotating shaft. Indeed, the lobe surface cooled more rapidly than its interior [START_REF] Hon | Emplacement and Inflation of Pahoehoe Sheet Flows: Observations and Measurements of Active Lava Flows on Kilauea Volcano, Hawaii[END_REF], which may lead to friction around the rotating shaft and disturbed the measurements. Future instrument refinement is required to protect the shaft from the outer cooling crust.

The results presented here, including the integration of viscosity field measurement and detailed textural analyses, are a first step toward quantifying the behavior of lava whilst accounting for all its components (melt, crystals and bubbles). Indeed, a major advantage of field measurements compared with those made in the laboratory is their ability to measure the behavior of a three-phase mixture of lava in its natural and original state. One of the difficulties in laboratory experiments is that sample preparation (i.e., melting the lava in the laboratory) leads to degassing and loss of volatiles, making accurate measurement of lava containing bubbles challenging, if not impossible. Direct comparisons between the laboratory experiments and measurements made in the field are therefore difficult for now. Recreating field conditions in the laboratory (in particular volatile content, oxygen fugacity, cooling rate and shear rate) is an ongoing challenge and a goal for future experiments [START_REF] Kolzenburg | The Rheological Evolution of the 2014/2015 Eruption at Holuhraun, Central Iceland[END_REF]. Therefore, field measurement using a rotating viscometer offers a promising solution to quantifying three-phase lava viscosity relations at high temperature. Future field viscometry holds the potential to evaluate the effects of three-phase mixtures on rheology, in particular the role of bubbles, and bubble-bubble or bubble-crystal interactions, over a large range of strain rate.

For this, future field viscometers need to be more versatile so as to measure small torques at low applied strain rates. This will allow the complete flow curve to be recorded (Pinkerton et al. 1995b). Future instruments also need to be lighter so that one operator can carry out rapid measurements from a vent to a flow front, and in very remote places. Such measurements can be envisaged on well-channelized narrow flows with stable levées or at breakouts at various distances from the vent (e.g., Belousov and Belousova 2018). With associated sample analyses and synchronous temperature measurement, field viscometry will be sufficiently precise to record the temperature-and time-dependent rheological transformations of lava and provide a benchmark for calibrating lava flow emplacement models. Measurements on various lava compositions, temperatures, and/or crystal and bubble content will also lead to improved quantification of lava rheology.

Conclusions

Viscosity measurements performed in the field by introducing a rotational viscometer directly into flowing lava, along with simultaneous sampling to allow the associated textural and chemical characteristics to be defined, are non-existent. We successfully measured the viscosity of several lobes that had an average of 380 Pa s for strain-rates between 1.6 and 5 s -1 and a temperature of 1144 °C. Although the full rheological curves could not be measured due to instrumental limitations and challenging field conditions, this result is in agreement with previous studies and confirms the range of apparent viscosities at unit strain rate for Hawaiian lavas. In contrast to previous studies, our measurements were made alongside simultaneous collection of quenched lava samples, which we thoroughly characterized by quantifying the chemical and textural characteristics. This provides a unique benchmark allowing field-based, laboratory-based and theoretical models to be compared and their accuracy, error and uncertainty to, for the first time, be assessed. Comparison of our measurements with viscosity estimated from a petrologic approach or from previous field and laboratory studies shows that the petrologic approach, when carried out with rigorous textural analyses of the collected samples, provides viscosity estimates that are in approximate agreement with field measurements. Quantification of the effect of crystals and bubbles is, however, extremely important if we are to correctly apply a model. Isothermal laboratory measurements prove to be the most precise method in describing the full flow curve, but cannot be applied directly to natural flows as crystallization takes place at higher temperatures and they cannot account for the effect of bubbles. We thus find that dynamic cooling experiments in conjunction with field measurements are needed to fully understand the lava emplacement rheology. Field viscometry accompanied by thorough analyses of quenched samples is the most promising approach to quantify evolving lava rheology in time and space. Collecting several measurements down a single flow system or at several flow units with different characteristics (temperature of effusion, cooling rate, composition, crystal and bubble content) will lead to proper benchmarking for lava flow modeling. Pinkerton et al. (1995) which are treated as a power law rheology. The grey curve represents the best power law fit and the black line is the linear Newtonian fit; the dashed lines represent the viscosity values as calculated from the petrologic approach (PA) using the chemical and textural characteristics averaged for all samples, and for samples quenched after viscosity measurements (R9, R10 ad R12). 3) FTIR method to determine the water content in the glass: Total water content was determined by Fourier Transform Infrared Spectroscopy, using a Bruker Vertex 70 spectrometer coupled with a Hyperion microscope system, housed at LMV. Spectra were recorded using a Globar light source, a KBr beamsplitter and a MCT (Mercury-Cadmium-221 Tellurium alloy) detector. Double polished samples were placed on a CaF2 window. Absorbance and background spectra of the glasses were obtained by acquisition of 300 and 100 scans, respectively. Beam size was about × 50 μm and spectral resolution was 4 cm -1 . Water concentration (OH groups + molecular H2O) was determined from the height of the absorbance band at ~3550 cm-1, using the Beer-Lambert law. We used the molar absorption coefficient provided by [START_REF] Mercier | Spectroscopic Analysis (FTIR, Raman) of Water in Mafic and Intermediate Glasses and Glass Inclusions[END_REF], 62.8 ± 0.8 L mol -1 cm -1 . Peak integration was performed using OPUS software.

Figure 1 :

 1 Figure 1: a) Schematic of the viscometer. b) Field viscosity measurements (Run 1) of the 61 G lava break out on November 23, 2016 (N19.349537, W155.048793).The two people most exposed to the hot lava are equipped with heat suits and goggles. One is holding the viscometer in position and the other is taking temperature measurements with an 80 cm long, K-type thermocouple with a 15 mm thick sheath.

Figure 2 :

 2 Figure 2: Results of testing the rotational field viscometer in the laboratory, using the viscosity standard N190000 from Cannon Instrument Company® at various temperatures; a) Example of raw data acquired in the laboratory; b) Flow curves obtained using both vanes; c) measured viscosity versus theoretical values of the standard viscosity liquid calculated using the relationship obtained from certified values at the experimental temperature. Horizontal error bars represent the viscosity variation across the plastic container due to the temperature gradient (if not visible they are smaller than the symbol). Vertical error bars represent the range of estimated viscosity over the different strain rates and instrument accuracy. The solid line gives the 1:1 relationship and dashed lines delimit an interval of ±50 Pa s. Squares are for measurements using the 80 mm long and 50 mm diameter shear vane (S80/25) and triangles are for measurements using the 60 mm long and 40 mm diameter shear vane (S60/20).

Figure 3 :

 3 Figure 3: In-situ measurements obtained at the 61G lava flow, on November 23, 2016 (N19.349537, W155.048793): a) Example of raw data acquired where horizontal double arrows show the range of constant speed over which torque was averaged; b) each data point is an average value over stable readings and error bars represent the fluctuation of the torque measurements (see supplement). Colored symbols are for individual runs treated as Newtonian, black symbols represent runs treated as having a power law rheology, and white symbols are the measurements fromPinkerton et al. (1995) which are treated as a power law rheology. The grey curve represents the best power law fit and the black line is the linear Newtonian fit; the dashed lines represent the viscosity values as calculated from the petrologic approach (PA) using the chemical and textural characteristics averaged for all samples, and for samples quenched after viscosity measurements (R9, R10 ad R12).

Figure 4 :

 4 Figure 4: Collected samples, from left to right field setting, hand sample, binary image of thin section scans (white is the glass, black are the vesicles and grey are the phenocrysts), SEM image at ×25 magnification: a) control sample (CN), collected on a breakout pāhoehoe lobe a few weeks old; b) example of a sample (C1) collected using a stainless steel tube inserted into the front of a pāhoehoe lobe and quenched rapidly in water, c) example of a sample (R12) collected from the lava attached to the vane after viscosity measurement and quenched in water. Red boxes indicate the thin section locations. Abbreviations are: Gl: Glass; Ol: Olivine, Pl: Plagioclase; Vs: Vesicles; LVs: Large vesicles formed during sampling (observed in quenched samples but absente from the control sample).

Figure 5 :

 5 Figure 5: The cumulative vesicle number density plot considers the vesicle number density per volume in mm -3 (Nv) with diameter greater than L (the equivalent diameter in mm) for all samples. C1a and C1b are both from sample C1 but from thin sections made perpendicular to each other. The orange and red squares represent the large vesicles related to the injection of air during sampling for R10 and R12, respectively.

  where data were selected at constant speed to extract the avearge torque. Black arrows show when the torque strnogly increased due to cooling preventing good measurment.

Table 1 :

 1 Field viscometry results. 718

	Exp. number	Vane	Time (s)	Rotations per minute	Angular velocity (rad/s)	Torque (N m)	+/-	Shear stress* (Pa)	+/-	Strain rate** (s -1 )	Viscosity Ŧ (Pa s)
	Run 1	80/25	7	7.6	0.80	0.24	0.03	769	90	1.6
			7	11.0	1.15	0.32	0.05	1030	194	2.3
	Run 2	80/25	23	10.5	1.10	0.24	0.02	767	79	2.2
	Run 8	60/20	3	18.6	1.93	0.20	0.04	1332	204	3.9
			3	21.5	2.25	0.27	0.02	1810	135	4.5
			4	23.8	2.50	0.32	0.02	2132	154	5.0
	Run 9	60/20	3	13.8	1.44	0.11	0.03	755	173	2.9
			7	19.7	2.06	0.17	0.03	1152	179	4.1
	Run 10	60/20	6	17.2	1.80	0.21	0.03	1359	177	3.6
			4	19.9	2.08	0.25	0.02	1664	116	4.2
			4	22.6	2.37	0.28	0.02	1842	135	4.7
	Run 11	60/20	6	12.3	1.29	0.16	0.04	1062	119	2.6
	Run 12	80/25	6	11.9	1.25	0.30	0.05	955	153	2.5
			11	16.3	1.49	0.40	0.03	1284	104	3.0
			5	17.6	1.84	0.43	0.04	1380	128	3.7
			5	23.0	2.40	0.55	0.06	1756	200	4.8
	*calculated with Eq. 1

**calculated with Eq. 3 and n = 1 Ŧ calculated with Eq. 4

Table 2 :

 2 Chemical composition of the 61G lava flow (bulk rock from sample R12) and average of the interstitial glass composition (average of 10 to 15 analyses per sample; numbers in parentheses are the standard deviations).

		Bulk rock (wt.%)	Glass (wt.%)
	SiO 2	50.82	51.95 (0.09)
	TiO2	2.45	2.77 (0.03)
	Al 2 O 3	13.71	13.33 (0.05)
	FeO *	11.18	11.84 (0.18)
	MnO	0.18	0.18 (0.01)
	MgO	6.91	6.27 (0.1)
	CaO	11.41	10.73 (0.08)
	Na 2 O	2.43	2.42 (0.03)
	K 2 O	0.62	0.50 (0.01)
	P 2 O 5	0.24	0.03 (0.01)
	Cr 2 O 3	-	0.02 (0.01)
	H 2 O tot	-	0.077(0.015)
	LOI	-0.68	-
	Total	99.36	100
	Mg #	0.38	0.35
	* total iron is presented as FeO
	Bulk rock oxide composition is normalized to
	a total of 100% but measured total is also
	reported
	LOI: Loss of Ignition
	H 2 O tot measure by FTIR see supplementary
	material

Table 3 :

 3 Textural analyses of the samples. Bulk density is the density of clasts measured with the Geopyc; density-derived vesicularity is obtained from bulk density and DRE density of 2847 kg m -3 . Standard deviations are based on two-to-four measurements per samples and are given in parentheses. φ v = fraction of vesicles (2D vesicularity); Ol+Px vol. % = percentage of olivine + pyroxene; R ol+px = aspect ratio of olivine + pyroxene; Plg vol. % = percentage of plagioclase; R plg = aspect ratio of plagioclase; φ xtl = vesicle-free crystal (phenocrysts + microlites) fraction.

	Sample	Bulk density (kg.m -3 )	Density-derived vesicularity φ	φv	Ol+Px vol. %	Rol+px	Plg vol. %	Rplg	φxtl
	CN	1441 (0.03)	49.37 (0.96)	0.52	3.6	1.7	6.7	2.8	0.21
	C1a	1030 (0.17)	63.87 (6.15)	0.53	3	1.9	2.8	3	0.12
	C1b	960 (0.02)	66.18 (0.53)	0.49	2.7	1.7	3.2	3.3	0.12
	R9	1294 (0.18)	54.52 (6.22)	0.42	6.3	1.8	5.4	3.3	0.20
	R10	1096 (0.12)	61.49 (4.21)	0.57	2.3	2	3.4	3.2	0.13
	R12	1145 (0.04)	59.75 (1.27)	0.44	3.6	1.7	4.7	2.9	0.15

Table 4 :

 4 Rheological parameters estimated via the petrologic approach using the chemical and textural characteristics averaged for all samples, and for samples quenched after viscosity measurements (R9, R10 ad R12).

			All samples	R9	R10	R12
	Melt Viscosity
	VFT parameters a	A	-4.55
		B	5901.7
		C	582.2
	Melt viscosity (Pa s) at 1144 °C a		330
	Effect of crystal on viscosity
	Mean crystal fraction b φ xtl	0.16	0.19 0.14	0.13
	Average crystal aspect ratio	R	2.4	2.6	2.6	2.3
	Maximum packing c φ m	0.51	0.50 0.50	0.52
	Relative viscosity (effect of crystals) d		2.0	2.6	1.9	1.8
	Melt + crystal viscosity (Pa s)		699	849	632	590
	Effect of bubbles on viscosity
	Mean vesicle fraction φ v	0.5	0.42 0.57	0.44
	Mean vesicle raidus		0.34	0.39 0.39	0.31
	Relative viscosity (effect of bubbles) e		0.31	0.41 0.24	0.38
	Melt + crystal + bubble viscosity (Pa s)		220	345	155	225
	a: calculated using glass composition (Table 2) via Giordano et al. (2008)
	b: bubble-free mixture considering phenocrysts and microphenocrysts
	c: calculated for rough particles according to Eq. 49 in Mader et al. (2013)
	d: calculated via Maron-Pierce (1956)
	e: calculated via Llewellin and Manga (2005) for Ca > 1

Supplementary material 1 1) Plot showing the convergence of Eq. 2 to Eq. 3 at Ro =infinity

	2)
	R i S80/25
	R i S60/20
	R o

Results of test of the rotational field viscometer in the laboratory, using the standard viscosity N190000 from Cannon Instrument Company® at various temperatures

	Temperature (°C)	T dev (°C)	Vane	Time (s)	Rotation per minute	Angular velocity (rad/s)	Torque (N m)	+/-	Shear stress* (Pa)	+/-	Strain rate** (s -1 )	Viscosity (Pa s)	+/-
	21	0.6	60/20	33.5	3.1	0.33	0.07	0.02	488	138	0.7		191
				32.3	6.0	0.63	0.15	0.02	971	143	1.4		103
				12.25	12.7	1.33	0.32	0.02	2137	147	2.9		50
				8.15	18.9	1.98	0.47	0.03	3146	175	4.4		40
	22	0.6	60/20	47.8	2.7	0.29	0.05	0.02	353	106	0.6		168
				26.1	6.9	0.72	0.14	0.02	944	132	1.6		83
				24.8	10.2	1.07	0.24	0.02	1602	128	2.4		54
				11.6	21.1	2.21	0.44	0.03	2914	178	4.9		36
	25	0.6	60/20	44	8.6	0.90	0.13	0.02	852	121	2.0		61
				38	9.5	0.99	0.14	0.02	947	123	2.2		56
				38	13.0	1.37	0.20	0.02	1349	138	3.0		46
				35	18.6	1.95	0.29	0.02	1942	136	4.3		31
	34	0.6	60/20	21	9.9	1.03	0.10	0.03	637	209	2.3		92
				38	7.7	0.80	0.07	0.03	495	206	1.8		116
				9	10.8	1.13	0.10	0.03	694	196	2.5		79
				27.5	15.4	1.62	0.15	0.03	1020	205	3.6		57
				35.6	19.7	2.07	0.21	0.03	1399	217	4.6		47
	20	0.6	80/25	11.5	6.9	0.73	0.41	0.03	1310	170	1.7		100
				60	12.2	1.28	0.74	0.04	2351	236	3.0		79
				19.5	8.6	0.90	0.52	0.03	1663	175	2.1		83
				25.5	6.4	0.67	0.40	0.02	1269	162	1.6		103
				22	5.5	0.57	0.32	0.03	1014	185	1.3		138
				44	2.3	0.24	0.13	0.03	423	169	0.6		297
	27	0.6	80/25	8.1	6.7	0.70	0.22	0.02	699	163	1.6		99
				16.6	6.6	0.69	0.22	0.02	708	154	1.6		96
				31.7	4.7	0.49	0.17	0.02	529	140	1.1		122
				46.4	8.5	0.89	0.28	0.02	902	161	2.1		77
				30.2	14.8	1.55	0.49	0.03	1568	171	3.6		47
	30	0.6	80/25	35	4.7	0.50	0.11	0.02	350	129	1.2		111
				21	7.4	0.77	0.18	0.02	563	137	1.8		76
				28	9.8	1.03	0.24	0.02	768	154	2.4		64
				21	13.4	1.40	0.33	0.02	1042	144	3.3		44
	40	0.6	80/25	36	10.2	1.07	0.13	0.02	417	136	2.5		54
				45	15.8	1.66	0.20	0.02	634	141	3.9		36
				39	12.7	1.33	0.16	0.02	512	141	3.1		45
				12.5	9.3	0.97	0.11	0.02	360	122	2.3		53
	*calculated
	with Eq. 1
	**calculated with Eq. 2 n=1

Acknowledgements

This research was financed by the French Government Laboratory of Excellence initiative n°

ANR-10-LABX-0006, the Region Auvergne and the European Regional Development Fund. MOC acknowledges the Auvergne fellowship for full support and the ClerVolc program 6 for field support. This is Laboratory of Excellence ClerVolc contribution n° 296. The authors gratefully acknowledge Matthew Patrick, the Hawaiian Volcano Observatory and Hawaiian Volcanoes National Park where work was completed under National park permit HAVO-2016-SCI-0064. The pole technique of LMV, including Jean-Louis Fruquière, Cyrille Guillot, Christophe Constantin, Jean-Luc Devidal, Jean-Marc Henot and Federica Schiavi, is acknowledged for instrument modifications, sample preparation and sample analyses. Field measurements were performed with the help of Alejandra Gomez-Ulla, Guillem Pinto-Gisbert, Matthew Burgess and Matthew Patrick who are greatly acknowledged. Additionally, Antonio Caponi, Alejandra Gomez-Ulla, Marie-Anne Ancelin, Valentin Gueugneau and Romina Severino are acknowledged for lab assistance during instrument calibration. Finally, Tamsin Mather (editor), Stefan Kolzenburg and Alessandro Vona (reviewers) are gratefully acknowledged for their thorough reviews and suggestions, which significantly improved the message of this contribution. References Applegarth, L J, H Tuffen, M R James, and H Pinkerton. 2013. "Degassing-Driven Crystallisation in Basalts." Earth-Science Reviews 116: 1-16. Ball, M., H. Pinkerton, and A. J L Harris. 2008. "Surface Cooling, Advection and the Development of Different Surface Textures on Active Lavas on Kilauea, Hawai'i." Journal of Volcanology and Geothermal Research 173 (1-2): 148-56. doi:10.1016/j.jvolgeores.2008.01.004. Barnes, H A. 1989. An Introduction to Rheology. Edited by Barnes H.A., Hutton J.F., and Walters K. Rheology Series 3. U.S. and Canada, Elsevier Science. Belousov, A., and M. Belousova. 2018. "Dynamics and Viscosity of 'A 'ā and Pāhoehoe Lava Flows of the 2012-2013 Eruption of Tolbachik Volcano, Kamchatka (Russia)." Bulletin of Volcanology 80 (6). doi:doi.org/10.1007/s00445-017-1180-2. Castruccio, A, A C Rust, and R S J Sparks. 2014. "Assessing Lava Flow Evolution from Post-