Magdalena Oryaëlle Chevrel

Andrew J L J L Harris

Alexian Ajas

Jonas Biren

Lucia Gurioli

Laura Calabro

Investigating physical and thermal interactions between lava and trees: the case of Kīlauea's July 1974 flow

Keywords: Lava-tree, lava channel, cooling rate, viscosity, pāhoehoe-ʻaʻā transition

come   L'archive ouverte pluridisciplinaire

Introduction

Lava flow emplacement dynamics are controlled by extrinsic parameters (topography, nature of the substrate, effusion rate, total volume emitted) and intrinsic properties (composition, temperature, vesicularity, cooling rate) which together control the viscosity and velocity of the lava (e.g., [START_REF] Chevrel | Lava flow rheology: a comparison of morphological and petrological methods[END_REF][START_REF] Harris | FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel[END_REF][START_REF] Kolzenburg | The rheological evolution of the 2014/2015 eruption at Holuhraun, Central Iceland[END_REF][START_REF] Rumpf | The influence of topographic roughness on lava flow emplacement[END_REF]. Under cooling-limited conditions [START_REF] Pinkerton | Methods of determining the rheolog-ical properties of magmas at sub-liquidus temperatures[END_REF][START_REF] Wilson | Mars review and analysis of volcanic eruption theory and relationships to observed landforms[END_REF], the run-out distance of a lava flow (i.e., the maximum length that can be reached at a given effusion rate) has typically been related to the rate of effusion (Walker 1973), to the total volume of material erupted [START_REF] Malin | Lengths of Hawaiian lava flows[END_REF] and to rheological changes due to the heat loss [START_REF] Pinkerton | Factor controlling the lengths of channel-fed lava flows[END_REF]. The evolution of the surface morphology of a channel-fed lava flow system commonly changes from a pāhoehoe-dominated proximal zone, to a medial transition zone with the formation of lava channels, to a distal zone of ʻaʻā [START_REF] Lipman | Aa flow dynamics, Mauna Loa 1984[END_REF]. The transition between pāhoehoe and a'ā has been the focus of many studies and is thought to be caused by an increase in strain rate or by a change in lava rheological properties including higher viscosity, development of yield strength and disruption of the cooled surface (e.g., [START_REF] Peterson | Transition of basaltic lava from pahoehoe to 'a'a, Kilauea volcano, Hawaii: field observations and key factors[END_REF][START_REF] Cashman | Cooling and crystalli-zation of lava in open channels, and the transition of pahoehoe lava to 'a'a[END_REF][START_REF] Sehlke | Pahoehoe to 'a'a transition of Hawaiian lavas: an experimental study[END_REF]. To characterise the changing properties that control lava down-flow behaviour, several techniques have been used. If the flow is active, lava properties may be directly measured in the field [START_REF] Lipman | Aa flow dynamics, Mauna Loa 1984[END_REF][START_REF] Moore | Preliminary estimates of the rheological properties of 1984 Mauna Loa Lava[END_REF]Crisp et al. 1994;[START_REF] James | Image-based measurement of flux variation in distal regions of active lava flows[END_REF][START_REF] Belousov | Dynamics and viscosity of 'a'ā and pahoehoe lava flows of the 2012-2013 eruption of Tolbachik volca-no, Kamchatka (Russia)[END_REF]. In the case of solidified flows, several studies have analysed the textural and chemical evolution of samples collected between the vent and the flow front to obtain down-flow thermal, textural and rheological evolution (e.g., [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[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][START_REF] Robert | Textural and rheological evolution of basalt flowing down a lava channel[END_REF][START_REF] Castruccio | The influence of effusion rate and rheology on lava flow dynamics and morphology: a case study from the 1971 and 1988-1990 eruptions at Villarrica and Lonquimay vol-canoes, Southern Andes of Chile[END_REF]. Other studies have modelled lava flows as function of cooling either with analogue materials (e.g., [START_REF] Hulme | The interpretation of lava flow morphology[END_REF][START_REF] Sakimoto | Channeled flow: analytic solutions, laboratory experiments, and applications to lava flows[END_REF][START_REF] Garel | An analogue study of the influence of solidification on the advance and surface thermal sig-nature of lava flows[END_REF] or via numerical models (e.g., [START_REF] Crisp | Influence of crystallization and entrainment of cooler material on the emplacement of basaltic 'a'a lava flows[END_REF][START_REF] Dragoni | The effect of crystallization on the rheol-ogy and dynamics of lava flows[END_REF][START_REF] Harris | FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel[END_REF][START_REF] Bernabeu | Modelling lava flow advance using a shallow-depth approximation for three-dimensional cooling of viscoplastic flows[END_REF][START_REF] Kelfoun | VolcFlow capabilities and potential devel-opment for the simulation of lava flows[END_REF].

The lava flow emplacement is also highly dependent on ground topography, substrate roughness [START_REF] Rumpf | The influence of topographic roughness on lava flow emplacement[END_REF]) and encountered obstacles [START_REF] Dietterich | Channel networks within lava flows: formation, evolution, and implications for flow behavior[END_REF]. [START_REF] Dietterich | Diverting lava flows in the lab[END_REF] examined how lava, using small basaltic flows created in the laboratory-effectively pouring molten basalt onto a synthetic topography-can be diverted by obstacles. Using the same method, [START_REF] Rumpf | The influence of topographic roughness on lava flow emplacement[END_REF] showed that the higher the ground roughness the lower the flow front velocity. Computational simulations were also used by [START_REF] Chirico | Lava flow hazard at Nyiragongo volcano, DRC 2 Hazard reduction in urban areas[END_REF] and [START_REF] Scifoni | Mitigation of lava flow invasion hazard through optimized barrier configuration aided by numerical simulation: the case of the 2001 Etna eruption[END_REF] to examine the effect of barriers on flow paths.

However, the dynamics of lava moving through forested areas remains rather unknown. Recently, [START_REF] Bernabeu | Laminar shallow viscoplastic fluid flowing through an array of vertical obstacles[END_REF] showed that a fluid moving through a grid of vertical cylindrical obstacles, as analogue for a lava flow inundating a dense forest, would form a wider and thicker flow than expected. Besides when the lava enter in contact with trees, energy will also be required to heat, dry and combust wood (e.g., [START_REF] Wagner | Calculations on forest fire spread by flame radiation[END_REF][START_REF] Babrauskas | Ignition of wood: a review of the state of the art[END_REF]. These studies suggest that the presence of cold, combustible obstacles in the flow may be added to the list of extrinsic parameters that potentially affect lava flows. Trees may thus play a role in both thermally and mechanically interacting with the flowing lava to exert a control on lava flow rheology, and thus also dynamics and morphology, as well as path and dimensions. However, although several studies exist that describe the morphology of casts left by trees due to the quenching effect on lava flowing around them, i.e., lava-trees (e.g., [START_REF] Finch | Lava tree casts and tree molds[END_REF][START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF][START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF][START_REF] Carveni | Lava trees and tree molds ( B cannon stones ^ ) of Mt[END_REF][START_REF] Parcheta | Hawaiian fissure foun-tains 1: decoding depositsepisode 1 of the 1969-1974 Mauna Ulu eruption[END_REF], no data are available on the effects of vegetation on lava flow propagation. There therefore remains an outstanding question: can trees cause enhanced lava cooling and mechanically interact with an advancing flow, thereby playing a role in changing flow dynamics and morphology (pāhoehoe versus 'a'ā) as compared with a forest-free scenario? To answer this question our study focuses on the July 19-22, 1974 lava flow of Kīlauea (Hawaii), a flow well-known for its abundance of standing lava-trees [START_REF] Hazlett | Geological field guide at Kilauea Volcano[END_REF]. We carried out detailed mapping and analysis of glassy samples collected from the inside and outside surfaces of the lava-tree casts and extracted flow temperature, plus crystal and bubble concentrations at each sampled site. These results were then used to estimate down-flow cooling and infer the corresponding increase in viscosity, which-for this case-appears no different to a tree-free emplacement case. Our analysis does, however, imply that some local mechanical interactions do occur with lava-trees, and allows us to propose a conceptual model and a viscosity threshold at which a tree may or may not turn into a lava-tree, and hence can or cannot interact. This point corresponds with, but does not cause, the pāhoehoe-'a'ā transition.

Lava-trees and tree-molds

It is well known that when low viscosity pāhoehoe invades a forest, often the lava does not topple the trees but rather cools rapidly against the cold, damp bark to form a cast of the tree [START_REF] Finch | Lava tree casts and tree molds[END_REF][START_REF] Jaggar | Volcanoes declare war: logistics and strategy of Pacific volcano science[END_REF][START_REF] Macdonald | Volcanoes and the sea -The geology of Hawaii[END_REF]. Following [START_REF] Finch | Lava tree casts and tree molds[END_REF], a tree inundated by lava will form a quenched cast of lava around it. The tree will then heat to combustion and burn slowly as the lava continues to flow around it [START_REF] Honda | Physico-chemical explanation for remeleting process of inner surface wall of Tainai tree molds located on the flank of Mt. Fuji[END_REF]. If the lava level remains constant then, once the tree has burnt away, a "hollow cylinder" [START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF] or "hole" [START_REF] Finch | Lava tree casts and tree molds[END_REF]) will be left in the lava flow with the same size and form as the tree root, trunk and branch system. The texture of the bark is often imprinted onto the interior of the cast surrounding the hole, so this is termed a "tree-mold" [START_REF] Finch | Lava tree casts and tree molds[END_REF][START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF]. If the lava level subsequently drops, then the tree cast will instead be left standing as a column (Fig. 1); this is a "lava-tree" [START_REF] Finch | Lava tree casts and tree molds[END_REF][START_REF] Macdonald | Volcanoes and the sea -The geology of Hawaii[END_REF][START_REF] Lockwood | Volcanoes global perspectives[END_REF]. In this case, rapid lowering of the lava level may cause the crust to scrape across the plastic outer surface of the cast to cause striations and tension-produced gash fractures [START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF]. In some cases, the lava may entirely overcome the height of the tree to form a closed lava-tree [START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF]. Where the lava did not reach the top of the tree, the hole left by the consumed trunk is exposed forming an open lava-tree. In such a case, one can have access to pre-lava inundation ground surface through the shaft (Fig. 1). In some places, tree trunks inside the cast may not burn completely, indicating that the lava inundation and drainage must have taken place in a short time [START_REF] Fiske | The December 1965 eruption of Kilauea Volcano, Hawaii[END_REF]; in others, the tree may be toppled and ingested horizontally into the lava or its crusts (basal or surface) to form horizontal, tubular molds [START_REF] Macdonald | Volcanoes and the sea -The geology of Hawaii[END_REF][START_REF] Carveni | Lava trees and tree molds ( B cannon stones ^ ) of Mt[END_REF].

Lava-trees have been used in various ways, such as to estimate maximum thicknesses reached by a lava flow [START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF][START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF], to retrieve flow direction [START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF], or to reconstruct the pre-eruption topography of the ground surface [START_REF] Jones | Proximal lava drainage controls on basaltic fissure eruption dynamics[END_REF][START_REF] Parcheta | Hawaiian fissure foun-tains 1: decoding depositsepisode 1 of the 1969-1974 Mauna Ulu eruption[END_REF]. Lava-trees have also been described as tools for "public educational purposes" [START_REF] Carveni | Lava trees and tree molds ( B cannon stones ^ ) of Mt[END_REF]) and, given their narrow cave-like properties, been the focus of speleological studies [START_REF] Bella | Tree mould caves within the framework of cave genetic classification[END_REF]. Additionally, charcoal left beneath a flow or in the tree molds [START_REF] Searle | A note on the formation of native iron and other effects associated with contact of basalt and carbonized wood at Auckland, New Zealand[END_REF]) can be analysed by C 14 for dating lava flow emplacement [START_REF] Kuntz | Radiocarbon studies of latest Pleistocene and Holocene lava flows of the Snake River Plain, Idaho: data, lessons, interpretations[END_REF].

All of these studies have described the tree-lava contact and the formation or structure of tree molds and lava-trees. However, no in-depth study exists on under which conditions lava-trees form and whether and how a forest might influence the emplacement dynamics and final form of a lava flow. [START_REF] Guest | The evolution of flow fields: observations of the 1981 and 1983 eruptions of Mount Etna, Sicily[END_REF], [START_REF] Harris | Effusive crises at Piton de la Fournaise 2014-2015: a review of a multi-national response model[END_REF] and [START_REF] Bernabeu | Laminar shallow viscoplastic fluid flowing through an array of vertical obstacles[END_REF] have suggested that trees can act to, or at least participate in, slowing down the advance of a lava flow. In addition, [START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF] have shown that if two trees are sufficiently close, a chilled crust may form between the trunks to form a wall or barrier of solid lava The July 19, 1974 lava flow Around 12:30 (Hawaiian Standard Time) on July 19, 1974, lava fountaining began from a small fissure at the base of the south wall of Keanakāko'i pit crater; with a second fissure opening a few minutes later north of Keanakāko'i [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF]. Shortly thereafter, two en-echelon fissures appeared on the south-eastern sector of the caldera floor, feeding a pāhoehoe sheet flow eastward and northward across the caldera floor [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF].

Around the same time, the fissures extended through the Ōhi'a forest to the west of Keanakāko'i to feed "a fast-moving fluid flow consisting of slabby pāhoehoe and 'a'ā that travelled quickly to the south and southeast" for a distance of 2 km [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF].

By 13:00, lava from these fissures began to cascade into Keanakāko'i to the north and Lua Manu craters to the east, forming a second, eastern, lava flow unit and cutting the chain of craters road (Fig. 2). Although the eruption continued until July 22, activity along the fissures to the east of Keanakāko'i ceased by 16:15 on July 19 [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF], having lasted around three and a half hours. [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] gives a volume for the south-eastern unit of 3.5 × 10 6 m 3 , for a mean output rate (given a 3-5 hour eruption duration) of 195-325 m 3 /s (280 m 3 /s for a 3.5 hour duration). It is the southeast lava flow of this fissure system on which we focus (Fig. 2). This south-eastern lava flow formed a "field of lava-trees" [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF][START_REF] Glatzer | Spectacular show at summit[END_REF] reports forest fires ignited by the passage of the lava, although these were "well under control by nightfall" on July 19. Instead, the "biggest headache" was "illegal parking (by sightseers) in the fire lanes", where "every possible place to park" was filled [START_REF] Glatzer | Spectacular show at summit[END_REF]. The damage caused to forest in this region due to the 1959 eruption of Kīlauea Iki and the 1974 eruption, as well as subsequent regeneration, has since been studied by [START_REF] Smathers | Hawai'i, the fires of life[END_REF]. The south-eastern lava flow is also the same flow as the one studied by [START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF] and [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF]. [START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF] used the lava-trees to approximate the maximum flow depth and thereby obtained maximum average flow velocities -1.56 m/s (near-vent) and 0.14 m/s (distally)assuming a viscosity of around 2000 Pa s. Velocities and viscosities are also available from features of super-elevation for a lava channel active during this eruption to the north of Keanakāko'i where lava "raced down a curving gully, with lava banking up against and spraying the outside of several bends" [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF]. In this channel, velocities were in excess of 8 m/s and viscosities were calculated by [START_REF] Heslop | Dynamics of a confined lava flow on Kilauea Volcano, Hawaii[END_REF] as being 80-140 Pa s.

In contrast, [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] examined the pāhoehoe-'a'ā transition in the south-eastern flow, which they estimated occurred at a distance of 1.5 km from the vent. In the present work, we examine whether the presence of trees affected downflow cooling, and hence rheological and dynamic changes in the channel-contained lava, as well as the relationship between the presence of lava-trees and the pāhoehoe-'a'ā transition.

Methods

Mapping and sampling

The lava flow area, surface morphology and lava-tree locations were mapped using digital elevation model (DEM) with spatial resolution of 10 m and a vertical accuracy of 1-2 m, as derived from TOPSAR (topographic synthetic aperture radar) by [START_REF] Pj | Quality of TOPSAR topographic data for volcanology studies at Kilauea Volcano, Hawaii: an assess-ment using airborne lidar data[END_REF]. The U.S. Geological Survey 7.5-minute (1:24,000-scale) topographic map sheet for the Kīlauea Crater quadrangle (1981 edition) was also used to support mapping, along with the cloud-free image satellite from GoogleEarth TM (image © 2018 DigitalGlobe -01/30/2016) in which individual lava-trees, as well as tree trunks lying on the flow surface, can be discerned. Locations of all lava-trees were added as a layer, and the density of lavatrees (number of lava-tree per unit area) was calculated from the number of objects within 200-m-long sectors along the length of the flow. The number of trees per unit area before the eruption was also estimated via counting treetops in areas of one hectare from an aerial image taken in 1965 (https://guides.library.manoa.hawaii.edu/aerials/digital -see online resource 1) and from an area aside of the flow using the GoogleEarth TM (image © 2018 DigitalGlobe -01/30/2016).

Additionally, field mapping and sampling were carried out in November 2016. During field work, 35 lava-trees within 200 m from the vent were located with hand-held GPS (WGS84) and their dimensions measured. Measured dimensions are marked on Fig. 1c and included: the diameter of the central hole, the width of the ring of solid lava surrounding the central hole, the height of the lava-tree above the current surface, and the depth of the central hole (i.e., distance from the highest flow level to the pre-existing ground surface). In addition, two sets of samples were collected. The first set corresponds to the sampling of nine lavatrees -one lava-tree every 200 m along a straight line in the centre of the channel between the vent and the flow front (T1-T9, Fig. 2 and3, Table 1), plus one lava sample located on the eruptive fissure itself (vent sample "V" in Fig. 2). This set of samples is used to study the lava chemical and textural variation down flow. The sampling involved the whole width of the cast at the top of the lava-tree, so that the hand sample included the cast surface that was in contact with the tree (hereafter called "tree-side") and the exterior of the cast that was in contact with the lava when it drained away (hereafter called "lava-side"). All dimensional measurements of Fig. 1c were also taken. Two to three thin sections were made from each lava-tree: one from the tree-side surface (hereafter labelled "-T"), a second from the lava-side surface (hereafter labelled "-L"), and-if possible-a third between the tree and lava sides (hereafter labelled "-M"), for a total of 26 thin sections (Fig. 3). The second set of samples was collected on the lava side of a single tree (T10, Fig. 4) to observe the lava chemical and textural evolution from the bottom to the top of the lava-tree cast. This lava-tree (T10) was selected as it was one of the tallest (1.92 m) found during the down-flow sampling, and showed six level mark-ers where lava had remained at a stable level for a short period of time during waning flow.

One sample was taken from each level (Fig. 4).

Chemical and textural analyses

The dense rock equivalent (DRE) density was determined from powdered sample fragments of known masses and using an Accupyc 1340 Helium Pycnometer that measured their volumes. The bulk density of the samples was measured using an envelope density analyser (Micromeritics Geopyc 1360), which measures the difference in the volume of a quasi-fluid medium with and without the sample embedded in that medium. As presented by [START_REF] Kawabata | Identifying multiple eruption phases from a com-pound tephra blanket: an example of the AD1256 Al-Madinah eruption, Saudi Arabia[END_REF] the Geopyc instrument can measure the envelope density of objects of different sizes and shapes. Prior to immersion, samples were dried in an oven for 24 hours at a temperature of 70 °C and then wrapped with thin laboratory parafilm®, to preserve external irregularities but to avoid the medium entering the porous sample. The standard deviation was estimated from five measurements to be ± 30 kg/m 3 by [START_REF] Thivet | Caractérisation magmatique du système superficiel du Piton de la Fournaise à travers l'étude des produits de l'éruption de Juillet[END_REF]. The density-derived vesicularity was then calculated using the relationship given by [START_REF] Houghton | A vesicularity index for pyroclastic deposits[END_REF].

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). The chemical composition of the glass and minerals were measured on polished thin sections via electronic microprobe CAMECA SX 100 (at 15 kV and a defocused beam of 20 µm at a current of 8 nA for glass and a focused beam at 15 nA for minerals). The temperature of lava during emplacement was obtained using the glass composition through applying the MgO-glass geothermometer of [START_REF] Helz | Geothermometry of Kilauea Iki lava lake, Hawaii[END_REF].

Backscattered electron (BSE) images were collected by scanning electron microscopy (SEM) at magnifications of 25× and 100× and converted into binary images that were used for extraction of vesicle size distribution (VSD). Following [START_REF] Shea | Textural studies of vesicles in volcanic rocks: an integrated methodology[END_REF], we considered only vesicles that were larger than 10 pixels, which corresponds to an equivalent diameter of 0.05 mm. Crystal content in the glassy area of the thin section at the contact with the surface of the tree side and lava side was estimated using two images per sample via ImageJ, considering only crystals larger than three pixels; that is, greater in length than 6 µm. The percentage of crystals was corrected for vesicularity and is therefore given for the vesicle-free mixture.

Estimation of lava viscosity and flow velocity

The viscosity of the three-phase lava mixture (fluid+crystals+bubbles) was estimated down the flow following the petrologic approach that considers the chemical, thermal and textural analyses of the samples (e.g., [START_REF] Pinkerton | Methods of determining the rheolog-ical properties of magmas at sub-liquidus temperatures[END_REF]Crisp et al. 1994;[START_REF] Guilbaud | Role of Syn-eruptive cooling and degassing on textures of lavas from the AD 1783-1784 Laki eruption, South Iceland[END_REF][START_REF] Robert | Textural and rheological evolution of basalt flowing down a lava channel[END_REF][START_REF] Chevrel | Lava flow rheology: a comparison of morphological and petrological methods[END_REF][START_REF] Chevrel | The viscosity of pahoehoe lava: in situ syn-eruptive measure-ments from Kilauea, Hawaii[END_REF]Rhety et al. 2017). This method was applicable only to samples containing sufficiently large (>1 mm 2 ) areas of glassy matrix where the composition of the quenched melt could be measured, temperature calculated and crystal content extracted. First the viscosity of the interstitial melt was calculated as a function of chemistry and temperature using the model of [START_REF] Giordano | Viscosity of magmatic liquids: a model[END_REF]. Second, the effect of the crystal cargo on viscosity was estimated following the method described in [START_REF] Mader | The rheology of two-phase magmas: a review and analysis[END_REF], which uses the equation of [START_REF] Maron | Application of Ree-Eyring generalized flow theory to suspensions of spherical particles[END_REF] and considering a crystal maximum packing calculated using the average aspect ratio of the crystals [START_REF] Mueller | The rheology of suspen-sions of solid particles[END_REF].

We thereby obtained the vesicle-free mixture viscosity. Third, the effect of bubbles on the mixture viscosity was estimated from the vesicle fraction (obtained from the vesicularity derived from density) via the equation given by [START_REF] Llewellin | Bubble suspension rheology and impli-cations for conduit flow[END_REF] for deformable bubbles. Using these viscosities, we calculated the instantaneous velocity of the lava with the Jeffreys equation [START_REF] Jeffreys | The flow of water in an inclined channel of rectangular section[END_REF]. For this, we considered the maximum and final depths of the flow as measured on each sampled lava-trees and an average underlying slope. Error propagation analyses on the estimation of lava viscosity and flow velocity are detailed in online resource 2.

Results

The July 19, 1974 lava flow and its lava-trees

We measured a length for the south-eastern lava flow of the July 19 flow field of 2200 m, with the unit covering an area of 60.55 × 10 4 m 2 in which a total of 598 lava-trees were mapped (Fig. 2). We note here that the flow generally contained lava-trees, and not treemolds, indicating a flow-wide lowering of the flow level as the eruption progressed.

Morphologically, the flow can be divided into five zones (Fig. 2 and5): near-vent pāhoehoe sheet flow, proximal channel, medial channel, distal channel and distal dispersed 'a'ā flow. Around the eruptive fissure, the lava is a smooth-surfaced "hummocky" [START_REF] Swanson | Pahoehoe flows from the 1969-1971 Mauna Ulu eruption, Kilauea volcano, Hawaii[END_REF]) pāhoehoe sheet flow over a surface with a slope of less than 0.5° (estimated from the lava flow surface). At 300 m from the fissure, multiple incipient channels arise within the pāhoehoe sheet flow. This is the zone of proximal channel which extends 200 m over a slope of 1.9° and feeds into the medial channel section at a distance of 500 m from the vent. Here, channels become well-formed across a zone of slightly steeper slopes (2.0°) and the flow becomes concentrated in two 20-30 m-wide main channels. As noted by [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF], here the two channels bifurcate around a topographic high and form streams of 'a'ā and thin, broken slabs of pāhoehoe surrounded by low initial levees of pāhoehoe; the channel margin is commonly only marked by a line of shearing between the static levées and the moving lava stream. At a distance of 900 m, the two channels coalesce into a single, broad (~50 m wide) and 500 m-long stream of 'a'ā and pāhoehoe slabs between broad initial levees of the same material. Down this distal channel section the slope is 0.75° and, while the occurrence of coherent pāhoehoe slabs within the channel decreases with distance, the percentage of surface area of 'a'ā increases (Fig. 2). After a distance of 1450 m, flow is of entirely 'a'ā. The zone of dispersed flow is 100 to 180 m wide, lies on a slope of 1.1° and feeds two 70-90 m-wide lobes with 2-3 m high flow fronts (Fig. 5).

Following [START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF], the thickness of the lava down flow was measured from the highest lava stand as given by the full lava-tree depth (i.e., value D in Fig. 1c). Then, following [START_REF] Jones | Proximal lava drainage controls on basaltic fissure eruption dynamics[END_REF] the final, lowest level was obtained from subtracting the lava-tree height above the current surface from the full lava-tree depth (i.e., value F = D -H in Fig. 1c). In the near-vent zone, lava-tree depths were in the range 1.1-7.7 m, with a mean of 3.2 m and a standard deviation of 1.3 m-this being the maximum flow level around the vent. Heights above the 1974 surface were between 0.42 and 1.95 m (1.1 ± 0.4 m), giving a final lava level of 0.5-6.5 m (2.1 ± 1.1 m). Down-flow, there was no systematic variation in lava level (Fig. 3). However, maximum levels (1.8-5.1 m, 1.9 ± 0.7, n. = 10), as well as minimum flow levels (0.2-2.4 m, 1.2 ± 0.7 m, n. = 10), were generally lower than around the vent. These values compare with the maximum flow thicknesses of [START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF] In terms of lava-tree density (Fig. 6), the proximal part of the lava flow contains the highest concentration (20-30 lava-trees/ha), with more than half of the lava-trees being located in the first 400 m near the vent (25% of the flow area) where pāhoehoe sheet flow dominated. As distance increases from the zone of pāhoehoe sheet flow, the density of lava-trees steadily decreases as the proportion of ʻaʻā to pāhoehoe increases (Fig. 6). There is a complete absence of lava-trees at the flow front, where 'a'ā is 100 %; and the last lava-tree is found 200 m from the front. At the flow front, trees in the process of being bulldozed can be found toppled onto the flow surface.

Morphology and dimensions of lava-trees

Near the vent and in the zone of pāhoehoe-sheet flow the lava-trees are, for the most part, well preserved. Most have conserved the initial tree trunk disposition of main trunk and branch systems and are lava-trees "en bouquet"; that is, they form a tightly grouped but upward splaying system of trunks and branches (Fig. 7a). Others are lone standing pillars either closed (Fig. 7b) or open (Fig. 7c) and are associated with the remnant (unburnt) trunks of the trees responsible for the casts (Fig. 7c). The heights of the lava-trees measured from the pre-flow surface to their top are between 1.1-7.7 m (in average 3.2 ± 1.2 m, n. = 42) and their casts have widths in the range 4-20 cm (10.12 ± 3.12 cm), with central hole (trunk) diameters of 6-70 cm (23 ± 11.5 cm). On the tree side of the cast the surface is smooth and tree bark imprints are well preserved. The lava side is rough, often with tension gashes, sharks teeth and striations that can be tied to descending slabs of pāhoehoe (Fig. 7d). All of these features indicate that the outer surface of the cast was still sufficiently plastic to "take the impressions of the crust" [START_REF] Moore | Lava tree molds of the September 1961 eruption, Kilauea volcano, Hawaii[END_REF] and to deform plastically [START_REF] Nichols | Superficial banding and shark's-touth projections in the cracks of basaltic lava[END_REF]. The descending pāhoehoe slabs remain in-place, surrounding the base of the lava-tree and uptilted towards the lava-tree and often forming a dome like structure around the lava-tree (Fig. 7e).

Lava accretion on the up-flow side of lava-trees and on the lava side of the casts was also often observed. The lava-tree of Fig. 4 is a good example. It is 1.9 m high, has a perimeter of 2.1 m, and is capped by a 20-cm-wide ring of rough-surfaced glassy pāhoehoe blebs, and has five level markers down the lava side of the cast. These markers are expressed as aligned trains of glassy globules and up-flow accretions. Glassy globules are stretched down flow (150° ESE) and plunge downwards at 22-25°. This is much steeper than the ground slope (2°), indicating rapid drainage in the down-flow direction. Up-flow accretions consist of larger blebs of rough surfaced pāhoehoe that wrap around the cast in the down flow direction, and have the same orientation and plunge as the trains of glassy globules (Fig. 4). This attests to maintenance of the high stand for a period of time to build the accreted cap, followed by a rapid decline in flow level with a series of brief stands at a single level to build the train-andaccretion sets. Given the 3-5 hour duration of the eruption, the total time for this process can be no more than an hour or two and with the transient declines and pauses being of the time scale of minutes.

In some places lava-trees are grouped (Fig. 7f) including two to four trees located at 1 to 2 meters of each other. In such a case, lava-tree casts merge to form a wall-like or barrierlike structure which can be aligned up or cross-flow (Fig. 7f). Their profiles are irregular with high points, marking the lava-trees, joined by lower ridges or "seams" [START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF] of accreted lava that join the lava-tree group together.

Textural and thermal analyses

The texture of the rock is porphyritic with phenocrysts of olivine within a glassy to microcrystalline matrix. Phenocryst content is 3 vol.% and does not change down flow. Microlites (<0.2 mm) include olivine and plagioclase, and their content increases from 0.6 vol.% near the vent to 33.1 vol.% near the front (Fig. 8 and9, Table 2). This gives a crystallization rate of around 15 vol.% per kilometre. As shown by the detailed sampling of the lava-side cast of T10, the microlite content appears to increase from 12 vol.% at the top of the lava-tree to 25 vol.% at the bottom (Table 3; Fig. 4).

Density and vesicularity derived from density (using a DRE density of 3030 kg/m 3 ) were in the range 1500-2000 kg/m 3 and 34-53 vol. %, respectively (Table 2, Fig. 10a). Only the sample from the tree side of T9 (near the vent) is markedly different with a density of 760 kg/m 3 and a vesicularity of 75 vol.%, representing at-vent gas-rich lava. Thereafter, there is no down-flow trend in vesicle content (Fig. 10a), but vesicularity of the lava next to the tree is always greater than that in the middle of the cast, attesting to a decrease in the bubble content of the erupted (and flowing) lava with time. The VSD analyses reveal that most samples contain one vesicle population between 0.05 and 5 mm, except T5-T that also comprises some vesicles larger that 10 mm and T1-T that includes two populations (Fig. 10b).

The chemical composition of the glassy matrix (online resource 3) decreases down flow from an Mg# of 36 at the vent to 28 near the flow front (at 1.9 km). For the tree side, this decrease in Mg# translates to cooling from 1139±2 °C at the vent and 1132±5 °C at 1.9 km.

For the lava side, we calculated cooling from 1140±2 °C to 1118±13 °C. The temperature decline estimated from the quenched lava against the tree is different to that for the lava side, where there is a cooling rate of 4 °C/km for initial lava quenched to the tree, and 10 °C/km for late-stage lava accreted to the outer portion of the cast (Fig. 9). These straddle the average cooling rate of 6 °C/km given for the same flow by [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] but show a lower cooling rate during initial flow than during late-stage flow. From the sampling of the lava-side cast of T10, we see that the temperature may change by 2 °C from level to level, but there was no systematic variation with depth during drainage (Table 3, Fig. 4).

Flow dynamics

The viscosity of the melt phase varies down flow as glass composition evolves and temperature decreases, so that the melt viscosity increases from 470 Pa s at the vent to 700 Pa s distally. Including the effect of the crystal fraction and deformable bubbles (Table 4), the viscosity of the three-phase mixture was estimated at 90-190 Pa s near the vent (considering the conditions at T9: 1140 °C, 3.6-13.5 vol.% crystals and 46-75 vol.% elongated bubbles). These viscosities are in-line with those of [START_REF] Heslop | Dynamics of a confined lava flow on Kilauea Volcano, Hawaii[END_REF] estimated from the super-elevation features in the near-vent channel active to the north of Keanakāko'i at the same time. Distally, the three-phase viscosity is up to 3600 Pa s (considering the conditions at T1: 1132 °C, 36 vol.% crystals and 50 vol.% bubbles). Our results give a viscosity (η in Pa s) with distance (Dist in m) trend of η = 160e 0.0017Dist with an R² of 0.94 (Fig. 11a).

This viscosity variation has a direct influence on the flow velocity (Fig. 11). At vent, at peak flow (maximum depth) velocities were as high as 9.8 ± 2.6 m/s (Table 4); similar to the maximum velocities obtained by [START_REF] Heslop | Dynamics of a confined lava flow on Kilauea Volcano, Hawaii[END_REF] but higher than what estimated by [START_REF] Moore | Estimates of lava-flow velocities using lava trees[END_REF]. With distance down flow, the increase in viscosity causes the velocity to fall to 0.3 ± 0.2 m/s at 1915 m. Declining output rate caused decreasing flow levels, so that the typical depth of late-stage flow near-vent (1.2 ± 0.7 m) gives a lower final velocity of 0.2-4 m/s at vent, waning to 0.04 m/s at 1915 m (Table 4). Given channel width of 30m, the maximum flow depth and the final flow level measured at T5, this translates to a decrease in effusion rate from 403 ± 250 m 3 /s to 15 ± 9 m 3 /s.

Discussion

Cooling effects of a forest

The 19 July 1974 lava flows advanced through a montane dry forest typical of Kīlauea's vegetation SW of Keanakāko'i crater [START_REF] Pratt | Terrestrial ecosystems[END_REF]. This zone receives less than 1000 mm of rain per year and experiences a dry season between mid-May and mid-October when rainfall is less than 100 mm/month [START_REF] Smathers | Hawai'i, the fires of life[END_REF]. The dominant forest tree is the Ōhi'a [START_REF] Pratt | Terrestrial ecosystems[END_REF], a pioneer species native to Hawaii that forms "the canopy for virtually all Hawaiian forests" [START_REF] Lincoln | Amy Greenwell Garden Ethnobotanical Guide to Native Hawaiian plants & Polynesian-introduced plants[END_REF]. All lava-trees identified in this study were of Ōhi'a. In the zone where the fissures feeding the south-eastern flow of the 19 July eruptive activity broke out, we estimated the number of trees before the eruption at 55 to 95 per hectare. Giving the average trunk radii of 10.5 cm and a typical depth of 3.2 m, this gives an area of 2.1 m² of cold bark and a volume of 11 m 3 of wood available to interact with the molten lava. The question is: are such tree densities, lava-tree contact areas and wood volumes sufficient to contribute a cooling effect to lava moving through the forest?

Lava flowing through the forest underwent an initial cooling rate of 4 °C/km which increased to 10 °C/km late in the eruption. This compares with cooling rates of 7 °C/km estimated for this same lava flow by [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] and for other lava flowing down channels in tree-free environments on Kīlauea and Mauna Loa (Crisp et al. 1994;[START_REF] Cashman | Cooling and crystalli-zation of lava in open channels, and the transition of pahoehoe lava to 'a'a[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] Robert | Textural and rheological evolution of basalt flowing down a lava channel[END_REF]. Thus, the cooling rates appear normal, so there was no discernible enhanced cooling effect due to the forest. The increase in cooling rate between the beginning and end of the eruption was likely due to the decrease in velocity in the channels with time as output rate feeding the channel system, and flow levels in the channels, declined.

Following [START_REF] Keszthelyi | Some physical requirements for the emplace-ment of long basaltic lava flows[END_REF], a lower velocity flow will undergo a higher rate of cooling per unit distance than a higher velocity flow. Simply, there is more time for lava moving at low velocity to cool as it passes from one location to the next than if it is moving at higher speed [START_REF] Harris | Effusion rate controls on lava flow length and the role of heat loss: a review[END_REF]. In our case, the increase in cooling rate is consistent with a decrease in velocity in the proximal section (Tables 3 and4).

The lack of a significant thermal interaction between the flowing lava and the trees, and hence no enhanced cooling effect, could be due to three reasons. First, the casts would have formed almost immediately and, as is the case for a tube roof [START_REF] Keszthelyi | Measurements of the cooling at the base of pahoehoe flows[END_REF], would have provided excellent thermal insulation to the lava flowing around the tree and thereby isolating the tree from lava flowing down the channel. Second, the velocity of the lava (2-10 m/s) and the small tree diameter (typically 0.2 m) would have meant that the flowing lava would have had a short interaction time (less than a tenth of a second) with the tree. Third, the forest had suffered damage due to fall-out during the Kīlauea Iki eruption of 1959 so that, as of 1962, the forest consisted of damaged Ōhi'a trees undergoing re-foliation, and lacked dense, mature undergrowth [START_REF] Smathers | Hawai'i, the fires of life[END_REF]. This, coupled with the fact that the eruption occurred during the middle of the dry season, would have meant that the moisture loads and vegetation densities would not have been as high as in a mature, undamaged forest. In the wet forest type that characterises the eastern edge of the flow and most of Kīlauea's east rift zone [START_REF] Pratt | Terrestrial ecosystems[END_REF], vegetation and water would have been more abundant still. Following [START_REF] Pratt | Terrestrial ecosystems[END_REF], such zones consist of a closed-canopy (i.e., high density) of Ōhi'a with an understory of dense tree fern. The dominant fern is the Ciboti-um species [START_REF] Pratt | Terrestrial ecosystems[END_REF], notably Hapu'u which can reach heights of 6-9 m, has "loose, fibrous bark on their trunks" [START_REF] Lincoln | Amy Greenwell Garden Ethnobotanical Guide to Native Hawaiian plants & Polynesian-introduced plants[END_REF]) and large (2-4 m long) fonds (Merlin 2005) which are loaded with water. Finally, locally, trees quenched the lava to solidify a zone of lava typically 10 cm wide. Considering the typical lava-tree height of 3.2 m, the quenched lava volume is 0.3 m 3 per tree, and around 180 m 3 within the flow. This volume is therefore rather small in comparison to the volume of the entire flow (3.5 × 10 6 m 3 ).

Mechanical interactions

The initial solidification of lava around the trees formed numerous lava-trees which remained in place as vertical pillars and provided a series of obstacles around which the lava had to pass. In the pāhoehoe sheet area of the flow, within 250 m of the fissure, about half of the trees have been changed into lava-trees, while further down flow, where the lava is of a'ā type, the lava destroyed all of the trees without creating any lava-trees.

Direction markers on the lava-side of the casts [START_REF] Lockwood | Lava trees and tree moulds as indica-tors of lava flow direction[END_REF] show that lava did have to move around each lava-tree. Three cases of diversion can be identified (Fig. 12): scattered (isolated), cross-flow grouped, and down-flow grouped. In the isolated case, there is an exceedingly local effect where flow paths diverge around the pillar-like obstacle which is just 0.06-0.6 m wide. In the case of the cross-flow grouped lava-trees, we have a barrier that can be up to 5 m wide (Fig. 7f). This impedes flow, causing flow to move laterally as well as to back up behind the barrier to build accretions that pile up on the up-flow side of the barrier. In the down-flow grouped case, a similar barrier to that formed in the cross-flow case is formed, but because of its down-flow alignment it has less effect on flow lines. However, the effect on diverging flow lines is greater than in the scattered case because enhanced accretion around the lava-tree group causes the obstacle to be wider (1-2 m) than in the case of an isolated lava-tree.

According to [START_REF] Dietterich | Diverting lava flows in the lab[END_REF], the velocity of the lava is halved after passing by an obstacle. It therefore seems that although we find the trees did not affect thermally the lava, the lava-trees, once formed, act as a network of solid pillars or obstacles around which the flow must move. Following [START_REF] Bernabeu | Laminar shallow viscoplastic fluid flowing through an array of vertical obstacles[END_REF] this can be viewed as a porous medium through which the fluid must move. The greater the number of lava-trees, the lower the porosity and hence the greater the effect on the dynamics for the fluid moving through the obstacle network. High obstacle densities cause the flow to be impeded and delayed, with the shape of the flow being changed (over the obstacle-free case), where fluid will spread and pile up in the porous zone [START_REF] Bernabeu | Laminar shallow viscoplastic fluid flowing through an array of vertical obstacles[END_REF]). Thus, a lava flow moving through a forest should be wider, thicker and shorter than for the same flow moving over a tree-free surface. This is also in agreement with previous observation of lava flows being thicker and slower where they branch or encounter obstacles [START_REF] Dietterich | Channel networks within lava flows: formation, evolution, and implications for flow behavior[END_REF]. This effect on lava flow advance could also be compared with the recent study of [START_REF] Rumpf | The influence of topographic roughness on lava flow emplacement[END_REF] where they find that flow front velocities are lower on rough terrain. However, they suggested that the slower advance rate is due to enhanced cooling through the higher surface contact area of an irregular substrate. Here we suggest that the thermal effect caused by the contact with the trees might be rather trivial and that the mechanical effect caused by vertical obstacles might be the main factor affecting the flow advances rate.

Interaction of lava with a forest might therefore affect lava flow hazard assessment because of the potential mechanical effect of the trees once turned into lava-trees. Following [START_REF] Bernabeu | Laminar shallow viscoplastic fluid flowing through an array of vertical obstacles[END_REF], further modelling is necessary to quantify the effect of vertical pillars as local lava diverter or flow delaying mechanism. As argued by [START_REF] Dietterich | Diverting lava flows in the lab[END_REF], "mitigation of lava-flow hazards must incorporate the dynamics of lava flow-obstacle interactions into barrier design". To go further, we suggest here that once and if the trees have been successfully transform into lava-trees, they may play the role of a "leaky" barrier [START_REF] Dietterich | Channel networks within lava flows: formation, evolution, and implications for flow behavior[END_REF]. Growing well-designed forests between potential source vents and vulnerable infrastructure may therefore offer a way to mitigate the risk caused by low viscosity lava flows. However, experiments should be carried out to test if and how a dense forest would affectively split the flow and slow down the advance. Scenarios including various types of trees and different tree arrangements need to be designed and tested.

Lava-trees and the pāhoehoe-'a'ā transition

There is a progressive down-flow transition from pāhoehoe to 'a'ā, beginning at a distance of 400-500 m from the vent, and is complete by a distance of 1500 m (Fig. 6). This is consistent with the findings of [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] who found the pāhoehoe-'a'ā transition occurring in this flow between 500 m and 1500 m. [START_REF] Soule | Examining flow em-placement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i[END_REF] ascribed the transition to crystallinity, with 'a'ā formation initiating when a critical crystallinity of 18 vol.% was reached. Likewise, we find that this transition initiates at 15-20 vol. %, which is the typical crystallinity between 400 m and 1000 m (Fig. 9), and a viscosity of 10 3 Pa s (Fig. 11).

We also find a coincidence between the incidence of lava-tree formation and the pāhoehoe-'a'ā transition which progressively declines from 20-30 lava-trees per hectare in pure pāhoehoe to zero in pure 'a'ā (Fig. 6). However, this is the pāhoehoe-'a'ā transition influencing the process of tree-formation; and not vice-versa. Simply, fluid pāhoehoe can flow between trees without toppling them, quenching to their trunks to allow lava-tree formation.

In contrast, 'a'ā bulldozes trees (as seen from the back-toppling cases at the flow front) knocking them down so that they do not remain as vertical structures around which the lava has to flow (Fig. 12). Thus, lava-tree formation (and flow impediment) is precluded. Instead, trees topple onto the flow to be rafted down flow to eventually be heated and burnt on the flow surface or rolled along in the basal clinker (Fig. 12). Thus, as lava progressively transitions to 'a'ā the incidence of lava-trees progressively declines. This decline follows a linear relationship whereby the number of lava-trees of per hectare (NT) is related to the percentage of lava as pāhoehoe (AA) through NT = -0.17(AA) + 17 (R² = 0.7). Considering that the number of pre-existing tree per hectare is constant, this means that in the pāhoehoe sheet area of the flow, within the 150 m from the fissure, nearly all the trees have been changed into lava-trees, while further down flow, where the lava is of a'ā type, all of the trees have been destroyed.

Re-creating the emission history from a single lava-tree Using the level markers on lava-tree T10 we can recreate the emplacement history of the lava flow. Initial flow entered the channel (17.5 m wide) in which T10 resided at an effusion rate of 305 m 3 /s (Table 3). This is 75 % of the value obtained for the single master channel at 920 m and into which this channel feeds (Table 4). Waning levels of flow in the channel reflected waning velocity and effusion rate, which dropped to a final value of 7 m 3 /s. This is now around 50 % of the value obtained for the master channel at the end of the eruption at 920 m.

Thus, it appears, the channel branch in which T10 resided took between 50 and 75 % of the total lava flux.

The distribution of effusion rate (E r ) with lava depth (D) shown in Table 3 has an exponential decaying form (! ! = 1.4! !.!" ! , R² = 0.97). If we distribute these effusion rates through time to fit to the total erupted volume (3.5 × 10 6 m 3 ) and eruption duration (3.5 hours) this gives a function for the duration at which each level (#) was maintained of duration =3.5! !!"## . If we distribute the volume through time using this function we obtain the timings shown in Table 5. These timings suggest that peak flow was maintained the first half of the 3.5 hour-long eruptive period, with waning emission occurring over the second half. During waning, flow level (lava thickness) was maintained for an increasingly short period of time with shut down being abrupt. This is consistent with the observation by [START_REF] Lockwood | Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawai'i[END_REF] that "discharge rates at all fissures gradually began to decline" from around 2 hours after the eruption began, with all activity having ceased another one hour and 40 minutes later.

Conclusion

We studied the relationship between lava flow dynamics and lava-tree formation for the case of the July 19, 1974 lava flow of Kīlauea. We mapped the lava flow morphology (pāhoehoe vs. 'a'ā) and the distribution of lava-trees along the flow path. We found that lava-trees are abundant (>20 lava-trees/ha) in the pāhoehoe area but sparse (<5 lava-trees/ha) in the 'a'ā section of the flow. We conclude that lava-trees can be formed when the lava flow is fluid enough; in such a regime lava-trees can form in spite of high flow velocities (~10 m/s). However, when the lava becomes too viscous, a living tree will be overridden and no lava-tree can be formed. This threshold occurs at a crystallinity of 15-20 vol% and viscosity of 10 3 Pa s where the pāhoehoe-'a'ā transition takes place so that trees begin to be bulldozed by the 'a'ā flow rather than flooded by pāhoehoe lava. If lava-trees form, markers on their outer surface allow the eruption history to be recreated, which in this case involved 1.5 hours of high effusion rate flow (at ∼400 m 3 /s) followed by 1.5 hours of waning flow and abrupt shutdown.

In this case, no enhanced cooling effect due to the presence of the trees was recorded, but this may be a function of the rapidity of the inundation and the type of vegetation, which was in poor health at the time of inundation. Mechanically, though, the lava was locally diverted by the array of solid cylinders formed by the lava-trees. Forests may thus have an effect in delaying flow advance and may be able to locally divert a flow from the expected path.

Further study on the effect of trees depending on their type will permit definition as to whether tree stands of certain types could be used, once turned into lava-trees, as local lava diverter or flow delaying mechanism. Note, though, that if the lava is more viscous than a critical threshold, that is here given at 10 3 Pa s, lava-trees will not form. Further work in this direction is clearly needed to assess the potential for vegetation barriers as a mitigation tool. Lava-trees (black points) were located based on aerial photography, TOPSAR images and ground truthing. "Sampled lava-trees" indicates location were rock samples were taken (T1 to T10; Table 1); "V" is the location of a sample taken at the vent. 2004) is also reported for comparison (7 °C/km R 2 =0.71). The linear regression for the crystallinity indicates 16.2 vol.% / km (R 2 =0.9).

Figure captions

  of 2.0-5.2 m (3.8 ± 1.1 m, n. = 6), and minimum thicknesses of 1.1-1.7 m (1.3 ± 0.2 m).

Figure 1 :

 1 Figure 1: Lava-tree formation (modified from Lockwood and Hazlett, 2010). a) Pristine forest

Figure 2 :

 2 Figure 2: Map of the south-eastern July 1974 lava flow showing the studied unit and location

Figure 3 :

 3 Figure 3: Sampling methodology along the lava flow and thin section positions on lava-tree

Figure 4 :

 4 Figure 4: Sampling of a single lava-tree (T10). From left to right: photo of T10 with location

Figure 5 :

 5 Figure 5: Evolution of width (empty circles) and thickness (full circles) of the studied lava

Figure 6 :

 6 Figure 6: Number of lava-trees per hectare and percentage of lava type (pāhoehoe vs. 'a'ā)

Figure 7 :

 7 Figure 7: Examples of lava-tree morphologies: a) lava-trees "en bouquet" near the vent; b)

Figure 8 :

 8 Figure 8: SEM images of the matrix crystallization down flow.

Figure 9 :

 9 Figure 9: Down flow variations in the temperature obtained from the glassy matrix on the lava

Figure 10

 10 Figure 10: a) Down flow variation of the vesicularity derived from density; b) Cumulative vesicle number density plot considering the vesicle number density per volume in mm -3 (Nv)

Figure 11

 11 Figure 11: a) Down flow variation in the velocity (diamonds) and viscosity (crosses). The

Figure 12 :

 12 Figure 12: Conceptual model for the interaction of lava and trees depending on the lava type.

Figure

  Figure 1

Table captions Table 1 :

 captions1 Sample location and description

Table 2 :

 2 Textural analyses down the flow

Table 3 :

 3 Textural analyses and viscosity estimation along a single tree (T10)

Table 4 :

 4 Viscosity and velocity estimation along the flow

Table 5 :

 5 Reconstruction of the effusion rate history

Table 1

 1

		GPS coordinates Latitude (N) Longitude (W)	Distance from vent (m)	Height of the lava-tree (m)	Description
	Vent	19°23'57.8710"	155°15'29.6680"	0	--	Glassy spatter t vent
	T9	19°23'57.3504"	155°15'29.5245"	15	2.6	Near the vent, sheet pahoehoe
	T8	19°23'51.2160"	155°15'28.5489"	205	3.8	Sheet pahoehoe zone
	T7	19°23'46.4604"	155°15'29.3742"	355	2	Gas-rich zone, shelly pahoehoe
	T6	19°23'39.4332"	155°15'27.0157"	580	1.88	Shear zone in the channel, slabby pahoehoe
	T10	19°23'39.0408"	155°15'24.6118"	600	1.92	Near T6, edge of the channel
						Tree with three trunks, in
	T5	19°23'29.2380"	155°15'23.2431"	920	5	channel of slabby
						pahoehoe
	T4	19°23'23.8272"	155°15'16.5212"	1180	2.7	Slabby pahoehoe
	T3	19°23'14.9604"	155°15'10.4030"	1505	2.7	Slabby pahoehoe to 'a'a lava
	T2	19°23'08.8764"	155°15'08.1606"	1705	--	Slabby pahoehoe to 'a'a lava
						First lava tree found
	T1	19°23'04.0596"	155°15'03.2633"	1915	3.2	since the front,
						essentially 'a'a lava

Table 2

 2 Mg-thermometer from 10 glass analyses per sample d considering the crystallinity and calculated via Mader et al. 2014 for phi max =0.48 (r=3.6) b caluclated via Giordano et al. 2008 e considering 40+/-0.1% of deformable bubble via Llewelin and Manga 2005 c including phenocrysts and microlites f Jeffreys equation with slope of 2° and density of 1.82

	Level	Sample	Distance from top (m)	Flow Depth (m)	Temperature (°C )a	Melt viscosity (Pa s)b	Crystallinity (vol.%)c	Relative viscosity crystal effectd	Viscosity of three phase mixture (Pa s)e	Velocity (m/s)f	Effusion Rate (m3/s)
						±		±		±		±		±		±		±
	High-Stand	A	0	3.12	1137	4	481	27	11.9	2.0	1.8	0.6	363	124	5.6	1.9	304	58
	1	B	0.4	2.72	1137	2	460	22	16.9	1.7	2.4	0.5	467	99	3.3	0.7	157	27
	2	C	0.75	2.37	1136	2	475	14	18.0	0.3	2.6	0.1	520	24	2.2	0.1	93	8
	3	D	0.95	2.17	1136	2	498	28	19.0	1.6	2.7	0.5	583	101	1.7	0.3	64	12
	4	E	1.25	1.87	1138	3	469	19	17.7	1.0	2.5	0.3	502	61	1.4	0.2	47	6
	5	F	1.5	1.62	1135	3	505	29	19.5	2.9	2.8	0.8	610	185	0.9	0.3	25	6
	Low-Stand	G	1.75	1.20	1137	2	491 20	25.1	2.9	4.4	1.0	921	213	0.3	0.1	7	2
	a

Table 3

 3 [START_REF] Giordano | Viscosity of magmatic liquids: a model[END_REF] using the glass composition and temperature reported in Table 2 b considering the crystal content and vesicularity reported in Table 2 and assuming phi max =0.48 (r=3.6) and deformable bubble via Llewelin and Manga 2005 c Calculated from Jeffreys equation with a density of 1740 kg/m3

	Tree	Zone	Distance from vent (m)	Max flow Depth (m)	Final Flow Depth (m)	Slope (°)	Melt viscosity (Pa s)a	Viscosity of three phase mixture (Pa s)b	Max velocity (m/s)c	Min velocity (m/s)c
								±		±		±		±
	T9-T	1: sheet flow	15	2.6	0.9	0.3	471	54	193	23	1.0	0.1	0.1	0.0
	T9-L	1: sheet flow	15	2.6	0.9	0.3	478	41	91	28	2.2	0.7	0.3	0.1
	T8	1: sheet flow	205	3.8	2.4	1.9	441	113	277	74	9.8	2.6	3.9	1.1
	T7	2: proximal channel	355	2.0	0.2	1.9	523	61	338	97	2.2	0.6	0.01	0.0
	T6	3: Medial channel	580	1.9	0.7	2.0	--	--	--	--	--	--	--	--
	T10	3: Medial channel	600	3.1	1.7	2.0	483	164	538	263	3.6	1.8	1.1	0.5
	T5	3: Medial channel	920	5.1	1.7	0.8	678	236	735	457	2.6	1.6	0.3	0.2
	T4	4. Distal channel	1180	2.7	1.4	0.8	--	--	--	--	--	--	--	--
	T3	5. Dispersed flow	1505	2.7	0.6	0.8	--	--	--	--	--	--	--	--
	T2	5. Dispersed flow	1705	1.8	1.8	1.1	--	--	--	--	--	--	--	--
	T1	5. Dispersed flow	1915	3.2	1.2	1.1	702	145	3634	1975	0.3	0.2	0.04	0.02
	a caluclated via

Table 4

 4 Time (hours) Effusion Rate (m3/s) Duration (h)

				Volume
				Emplaced (m3)
	1.8	403	1.8	2.55 × 106
	2.6	209	0.9	6.57 × 105
	3.1	124	0.4	1.95 × 105
	3.3	85	0.2	6.68 × 104
	3.4	63	0.1	2.48 × 104
	3.4	34	0.05	6.64 ×103
	3.5	9	0.03	8.93 × 102
	Total	280	3.5	3.5 × 106

Table 5

 5

Acknowledgements

The authors gratefully acknowledge the support of 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 (Christophe Constantin, Jean-Luc Devidal, Jean-Marc Henot, Mhammed Benbakkar and Claire Fonquernie) is acknowledged for sample preparation and sam-ple analyses. Fieldwork was performed with the help of Scott Rowland and Alejandra Gomez-Ulla who are greatly acknowledged. Additionally, we thank Taeko Jane Takahashi at the HVO library for helping us track down the internal and newspaper reports for the events of 19 July 1974. Finally, T. Gregg, A. Soule (reviewers) and H. Dietterich (editor) are greatly acknowledged for their thorough reviews and comments, which improved the quality of this work.

Funding information

This research was financed by the Agence National de la Recherche through the project LAVA (Program: DS0902 2016; Project: ANR-16 CE39-0009, http://www.agence-nationale-recherche.fr/ Projet-ANR-16-CE39-0009). This is ANR-LAVA contribution no. 7. Fieldwork was supported by the Laboratory of Excellence ClerVolc pro-gram 6, contribution no. 316. MOC acknowledges the Auvergne fellow-ship for support.


James Gillespie
email: gillespie-j10@email.ulster.ac.uk

Iñaki Rañó

Nazmul Siddique
email: nh.siddique@ulster.ac.uk

José Santos
email: ja.santos@ulster.ac.uk

Mehdi Khamassi
email: mehdi.khamassi@isir.upmc.fr

Reinforcement Learning for Bio-Inspired Target Seeking

Keywords: Braitenberg Vehicles, Reinforcement Learning, Source seeking

Because animals are extremely effective at moving in their natural environments they represent an excellent model to implement robust robotic movement and navigation. Braitenberg vehicles are bioinspired models of animal navigation widely used in robotics. Tuning the parameters of these vehicles to generate appropriate behaviour can be challenging and time consuming. In this paper we present a Reinforcement Learning methodology to learn the sensori-motor connection of Braitenberg vehicle 3a, a biological model of source seeking. We present simulations of different stimuli and reward functions to illustrate the feasibility of this approach.

Introduction

Everyday experience shows that animals are extremely effective at moving in their natural environments which makes them an excellent model to implement robust robotic movement and navigation [START_REF] Floreano | Robotics and neuroscience[END_REF]. Bio-inspiration has played an important role in robotics for the design of locomotive systems [START_REF] Menciassi | Bio-inspired solutions for locomotion in the gastrointestinal tract: background and perspectives[END_REF], locomotion control [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF], and control of steering [START_REF] Webb | A Spiking Neuron Controller for Robot Phonotaxis[END_REF]. This paper focuses on learning animallike steering behaviour using Braitenberg vehicles, a well known model of insect tropotaxis [START_REF] Braitenberg | Vehicles. Experiments in synthetic psycology[END_REF], i.e. movement of insects towards or away from a stimulus [START_REF] Fraenkel | The orientation of animals. Kineses, taxes and compass reactions[END_REF]. While previous works used experimentally tuned models or relied on time consuming evolutionary strategies [START_REF] Bicho | The dynamic approach to autonomous robotics demonstrated on a low-level vehicle platform[END_REF][START_REF] Mondada | Evolution of neural control structures: some experiments on mobile robots[END_REF][START_REF] Salomon | Evolving and optimising braitenberg vehicles by means of evolution strategies[END_REF], this work investigates how robots can learn these tropotaxis controllers from experience using Reinforcement Learning techniques.

How animals navigate in the presence of stimuli has been a subject of research for several decades [START_REF] Fraenkel | The orientation of animals. Kineses, taxes and compass reactions[END_REF][START_REF] Gallistel | Navigation: Whence our sense of direction?[END_REF], yet the mechanisms generating these behaviours are not fully understood even for very simple animals. One natural mechanism to control movement is tropotaxis, which relies on sampling the stimulus in the environment at two symmetrically placed sensing organs (eyes, ears. . . ), and can be exploited to control robot movement. The values of the stimulus are used by the animal to steer or control its navigation, via some neural wiring between the sensing organs and the motor effectors. The concrete internal wiring defines the behaviour of the animal generating positive or negative taxis (movement towards or away from a stimulus), as captured in the seminal work of Valentino Braitenberg [START_REF] Braitenberg | Vehicles. Experiments in synthetic psycology[END_REF], where taxes are modelled through simple vehicles.

Braitenberg vehicles have been used in robotics for several decades from their most basic formulation [START_REF] Lilienthal | Experimental analysis of smelling braitenberg vehicles[END_REF], to extensions with fuzzy controllers [START_REF] Yang | A fuzzy-braitenberg navigation strategy for differential drive mobile robots[END_REF] or neural networks [START_REF] French | Introducing neuromodulation to a braitenberg vehicle[END_REF]. Therefore, they are widely used to implement bio-inspired robotic behaviours, especially when the motion relies on unconventional sensors, for example sensors which are not providing distance readings [START_REF] Lebastard | Underwater electronavigation in the dark[END_REF][START_REF] Lilienthal | Experimental analysis of smelling braitenberg vehicles[END_REF][START_REF] Salumäe | Against the flow: A braitenberg controller for a fish robot[END_REF]. Early research in the field aiming at enabling learning of animal-like movement relied on techniques like evolution strategies [START_REF] Salomon | Evolving and optimising braitenberg vehicles by means of evolution strategies[END_REF] and genetic algorithms [START_REF] Mondada | Evolution of neural control structures: some experiments on mobile robots[END_REF]. These works obtained neural architectures of Braitenberg vehicles to navigate environments with no collisions through an iterative approach to optimise a fitness functions of the distance to obstacles and forward speed. A big drawback of these strategies is the amount of time necessary for the algorithms to find the optimal weights of the neural network to control the velocity of the wheels. Moreover, having such a neural controller makes the analysis of the sensori-motor connection extremely difficult. Another early approach to bio-inspired robotic steering used the Dynamical Systems Approach to behaviour generation for obstacle avoidance combined with a Braitenberg vehicle for target seeking [START_REF] Bicho | The dynamic approach to autonomous robotics demonstrated on a low-level vehicle platform[END_REF]. Although the work presented is highly effective for robot navigation, the parameters of the control mechanism are selected experimentally, which means they might not generate the best possible behaviour.

Interestingly, recent works using Braitenberg vehicles have also developed models which do not use learning techniques to achieve their target seeking behaviour and can be sub-optimal. For instance, a basic implementation of Braitenberg vehicles to determine the speed of the wheels of a smelling robot [START_REF] Lilienthal | Experimental analysis of smelling braitenberg vehicles[END_REF] was used with a simple dynamical normalisation of the measured values of the sensors. All these works rely on experimental results, and while mathematical models of Braitenberg vehicles have been presented in [START_REF] Rañó | A model and formal analysis of braitenberg vehicles 2 and 3[END_REF][START_REF] Salumäe | Against the flow: A braitenberg controller for a fish robot[END_REF], it can be challenging to tune the internal parameters of the vehicle, i.e. the controller, as formal proofs of stability are missing. Moreover, an outstanding open question in Braitenberg vehicles is how to define the relationship between sensors and motors in an optimal way. The main contribution of this paper is presenting a reinforcement learning based methodology to perform learning in bio-inspired steering, as a way to answer the question of optimal Braitenberg vehicle design. Obtaining an optimal controller depends on the selected definition of optimality, but also on the stimulus the vehicle perceives. We present results for two common families of stimuli, namely following the inverse square law, and the inverse distance law. As we will see, known theoretical results of Braitenberg vehicles help analysing and interpreting the results from the reinforcement learning process. The rest of the paper is organised as follows. Section 2 presents a brief introduction to the model of Braitenberg vehicle 3a, and states the problem of tuning the sensorimotor connection as a reinforcement learning problem. Section 3 presents the experimental results obtained for the simulation of these vehicles using the types of stimuli mentioned and several reward functions defining optimality criteria. The paper ends presenting some conslusions and future work in section 4.

Reinforcement Learning in Braitenberg vehicles

Braitenberg vehicles are well known models of animal behaviour, used in robotics as a simple way of implementing avoidance and target seeking behaviours. The lack of a mathematical formalism for the vehicles did not hamper their usage but made their implementation a trial and error process, where intuition and experience played a big role. The development of a mathematical model of the closed-loop system of Braitenberg vehicle 3a [START_REF] Rañó | A steering taxis model and the qualitative analysis of its trajectories[END_REF] as a non-linear differential equation allowed using dynamical systems theory to analyse characteristics of the solution trajectories. Oscillatory [START_REF] Rañó | Results on the convergence of braitenberg vehicle 3a[END_REF] and unstable [START_REF] Rañó | A steering taxis model and the qualitative analysis of its trajectories[END_REF] trajectories can be found among the solutions of this non-linear dynamical system. Moreover, their formalisation and the model equations allow to get a deeper understanding of how these controllers work, and enable making better informed guesses on the parameters without the need for experimentation. However, as we already stated, one outstanding open question is finding under which conditions the behaviour of Braitenberg vehicle 3a is stable close to a given equilibrium point, since the linear stability test provides no information [START_REF] Rañó | Results on the convergence of braitenberg vehicle 3a[END_REF], while finding a Lyapunov function is a non-trivial problem. Therefore, our goal with this work is learning a stable non-linear controller that achieves target seeking for a Braitenberg vehicle 3a.

Figure 1 shows Braitenberg vehicle 3a in the proximity of a light source (although to derive its mathematical model the stimulus can be of any nature and does not need to come from a source). The wheels of the vehicle are used to abstract the locomotive systems of animals focusing on the steering level of motion [START_REF] Arechavaleta | An optimality principle governing human walking[END_REF]. In fact, similar steering models have been used to understand human motion. As shown in the figure, the vehicle has two sensors connected to the ipsilateral wheel (the same side wheel) in a decreasing (inhibitory) way, i.e. in the sense the higher the stimulus value, the slower the wheel turns. Intuitively the vehicle turns towards, and approaches, the light source, performing a 'hillclimbing' on the stimulus. The closer the vehicle gets to the source, or maximum, the lower its velocity, and it will eventually stop near the peak when the speed of both wheels is zero. Because intuitively the motion converges to the maximum, these vehicles were used in real robots and simulated artificial agents without mathematical formalisation of their behaviour. These works assumed the vehicle converged regardless of the stimulus provided, if the connection between the sensors and the wheels was appropriately tuned.

We will present the steps to model the vehicle shown in figure 1, with a wheelbase d and distance between the sensors δ. Assuming the vehicle is immersed in a stimulus that does not change over time, the stimulus itself can be modelled as a non negative function S(x) of the position in the environment or the ve- hicle workspace x = (x, y) ∈ W ⊆ R 2 . The stimulus is considered to simply be a position dependant scalar function, like light or sound intensity, for instance. Furthermore, we will assume S(x) is of type C ∞ in W, i.e. its value changes smoothly in the environment. If the stimulus is originated by a point-like source located at x S ∈ W the function S(x) will reach a maximum at this point as the stimulus intensity decays with distance, i.e. S(x S ) ≥ S(x)∀x ∈ W, which also implies ∇S(x S ) = 0 and y T D 2 S(x S )y < 0∀y ∈ R 2 because the function is smooth, i.e. the gradient of S(x) vanishes and the Hessian is a negative definite matrix at x S . Without loss of generality, as figure 1 shows, we can assume the source is located at the origin of a reference system, i.e. x S = 0, and the position of the vehicle is referred to that coordinate system.

According to the qualitative model proposed by Braitenberg, the connection between the sensor readings 's' and the left and right wheel velocities v L/R is direct (ipsilateral) and inhibitory. We will model this connection as a smooth function F (s), fulfilling the following criteria: (i) the motion of the vehicle is never backwards, i.e. F : + ∪ {0} → + ∪ {0}; (ii) the vehicle stops at the source, i.e. F (S(0)) = 0; and (iii) an inhibitory connection implies that the derivative of F (s) is negative for all stimulus values s, i.e. F (s) < 0. In the literature of Braitenberg vehicles for source seeking the function is selected to be an experimentally tuned linear function, or an Artificial Neural Network, however it can be any function fulfilling these conditions. Under the above conditions, the behaviour of the vehicle can be approximated as:

ẋ = F (S(x)) cos θ ẏ = F (S(x)) sin θ θ = - δ d ∇ x F (S(x)) • êp (1)
where x = (x, y) is the midpoint between the sensors, ∇ x F (S(x)) is the gradient of the function composition, '•' represents the dot product, and êp = [-sin θ, cos θ] is a unit vector perpendicular to the vehicle direction. While condition (ii) forces the dynamical system to have an equilibrium point, condition (iii) is related to the stability of the equilibrium point. Therefore the functions sought using Reinforcement Learning should fulfil these conditions for the vehicle to perform positive taxis, i.e. it should have negative slope and it should be zero when the vehicle is at the source.

Simulations of Vehicle stimulus and connection

The mathematical model introduced above is the tool to obtain analytic results on the behaviour of the vehicle, like stability analysis, features of the trajectories, et cetera. However, because its application in a reinforcement learning set-up would require computing the derivative of the function to learn F (s), we opted for a simpler statement of the learning problem. Hence, in our simulations we use the straightforward implementation of vehicle 3a, i.e. evaluating the stimulus at two points and computing the speeds for the left and right wheels. Given the position of the right and left sensors x r and x l , respectively, we can compute the velocities for the right and left wheel as v r = F (S(x r )) and v l = F (S(x l )), where F (s) is the sought function and S(x) is a selected stimulus function. The left and right sensor positions can be easily computed from the vehicle pose, i.e. the midpoint between the sensors, x and the orientation of the vehicle, θ, as x r = xδ 2 êp and x l = x + δ 2 êp , where the êp = [-sin θ, cos θ] is obtained from the vehicle orientation θ. The velocities of the wheels can be converted into forward speed v and turning rate ω of the vehicle as:

v = F (S(x r )) + F (S(x l )) 2 ω = F (S(x r )) -F (S(x l )) d ( 2
)
where it is worth remembering that the positions of the sensors depend on the vehicle pose, i.e. x r = x r (x, θ) and x l = x l (x, θ), and, therefore, the movement of the vehicle depends on its pose (x, θ) through the stimulus function and the sensori-motor connecting function F (s). Given these velocities, the trajectory of the vehicle can be obtained as the solution of the system of non-linear differential equations:

  ẋ ẏ θ   =   cos θ 0 sin θ 0 0 1   v(x, θ) ω(x, θ) (3)
which has to be integrated numerically, and in our experiments we used the Euler method with a fixed time step of 0.05. It is worth remembering that the connection function F (s) must be decreasing for the vehicle to be of type 3a, otherwise instead of target seeking the vehicle will move away from the source (ipsilateral increasing connection corresponds to vehicle 3b, negative taxis). However, as we will see in the next section this constraint cannot be imposed in the learning process, but good solutions are found using appropriate reward functions.

The reinforcement learning problem

The problem of learning the target seeking behaviour for Braitenberg vehicle 3a given a stimulus function can be casted into finding a non-increasing function of the stimulus to compute the velocities. Given the stimulus S(x) defined in some environment W ⊆ 2 and an initial pose of the vehicle (x 0 , θ 0 ) ∈ W × S 1 the trajectory followed depends on the connection function F (s). Since the trajectory unfolds in time (x(t), θ(t)) but also depends on the initial pose and F (s) we will write (x(t, x 0 , θ 0 , F ), θ(t, x 0 , θ 0 , F )). For each pose of the vehicle we can define a scalar reward function r(x, θ) measuring how good being at that state is, and because the trajectory depends on the initial pose and the connecting function, the reward along a given trajectory is therefore a function of the initial pose and the connecting function too, that is r(x, θ) = r(t, x 0 , θ 0 , F ). Theoretically, the function F (s) can be approximated using any type of function approximation, however, for this work we opted to use a Radial Basis Function neural network of the form:

F (s) = N i=1 w i φ i (s) ( 4
)
where N is the number of basis functions, w i is the weight of the i-th neuron, and the radial basis functions φ i (s) were Gaussian kernel functions centred at fixed equidistant positions s i within the range of the stimulus (from s = 0 to s = S(0)). If we denote the weight vector of the network by Φ = (w i ), the reinforcement learning problem can be stated as maximising the following total return:

R[Φ] = (x0,θ0)∈W×S 1 dx 0 dθ 0 ∞ 0 r(t, x 0 , θ 0 , Φ)dt (5)
where we need to integrate for all the initial conditions, i.e. the points in the workspace x 0 and all orientations in the unit circle θ 0 ∈ S 1 , but we also need to integrate over the whole trajectory. Because this total return function depends on the solution of a non-linear dynamical system it is impossible to evaluate these integrals. Moreover, integrating over the whole trajectory is obviously not feasible, so we defined a finite time t f to run the simulations, changing the upper integration limit in equation ( 5). This does not solve the problem of evaluating the return, however we can use the sampling trick to estimate through simulations the value of the integral by randomly sampling the space of poses of the vehicle. We can use roll-outs from random initial poses to estimate the gradient ∇ Φ R[Φ] and perform a hill climbing on the return to find the optimal weight vector of the RBF network using:

Φ k+1 = Φ k + α k ∇ Φ R[Φ k ] (6)
where a k = a0 k , a 0 is the initial learning rate, and, as already stated, the gradient was estimated using roll-outs through small random perturbations of the parameters Φ using central finite differences.

Although in our case the environment of the vehicle was the whole plane W = 2 , to simplify things further the domain of the integral in equation ( 5), i.e. the domain in which the initial conditions are selected for the sampling process, was defined to be a square region around the origin. Moreover, the initial angular directions of the vehicle were selected to be pointing towards the source within a ±90 • range, i.e. θ 0 ∈ [θ t -π/2, θ t + π/2], where θ t = arctan y0 x0 -π.

Experimental results

To learn the sensori-motor connection for a Braitenberg vehicle using Reinforcement Learning we defined and tested several reward functions. Because the goal of the vehicle is to reach the source of a stimulus, i.e. where it takes its highest value, an immediate candidate for the reward function would be the stimulus itself. Provided the vehicle has access to its pose, and knowing that the source is at the origin of the reference system, an alternative reward would be a function of the distance from the vehicle to the source and its relative heading. We simulated four different reward functions, two dependent on the stimulus and two on the pose of the vehicle, and to fulfil assumption (i) in section 2 the reward was given only when the movement of the vehicle was in the forward direction.

The reward functions used are:

1. The stimulus itself, and because the vehicle obtains readings from both sensors, we use as reward their sum, i.e. r = S L + S R . 2. The previous reward tries to maximise the stimulus value, but for a stimulus source it is important that the vehicle heads in the right direction. We defined a second stimulus based reward that accounts for the heading direction by trying to make the value in both sensors identical. The reward function accounts for the value but penalises directions perpendicular to the source r = S L +S R 1+(S L -S R ) 2 . 3. In the simulations we have access to the pose of the vehicle, (x, θ), which allows us to define additional reward functions. Using the pose this reward function consists of a linear combination of terms accounting for the proximity of the vehicle to the source and its heading. The selected function was r = a 1+α(θ,θ T ) 2 + 1 1+||x|| 2 , where a = 3 represents the relative importance of heading vs. the distance, ||x|| is the distance to the source, and α(θ, θ T ) is the angular distance in the range [-π; π] between the robot heading (θ) and the desired heading θ t = arctan y x -π. 4. Our early experiments with the reward functions defined above showed the vehicle's trajectories stopped way before reaching the source. We thought it was due to the limited/bounded rate of growth of these reward functions close to the source. Therefore, we included as a reward function r = 1 ||x|| 2 , but since it is singular at x = 0, we defined an threshold distance such that if ||x|| < , r = 1 2 . These four reward functions were used for the two general types of stimulus defined below. It is worth noting that in all the simulations the robot initial heading was randomly selected within ±π/2 radians in the direction of the stimulus.

Inverse-square law stimulus

Probably the best known example of Braitenberg vehicles is the one implementing phototaxis using light sensors. With a light source placed at some height h 0 above the origin of a reference systems of a ground plane, and since light intensity follows the inverse-square law, the stimulus as a function of the position x will be S(x) = In our first set of simulations we used the stimulus sum reward function, and the results are shown in Figure 2. We ran 10 learning experiments per reward function with random initialisations of an RBF network. Only results of the successful simulations are shown in the figures: Figure 2 depicts simulated trajectories on the left and the corresponding F (s) functions approximated by the RBF networks while Figure 3 shows an example of the evolution of the reward as a function of time. The successful simulations are those which fulfil the conditions stated above, i.e. taking positive values (the vehicle moves forward), having a negative slope (the behaviour is positive taxis), and they become zero close to the stimulus maximum (in our case S(0) = 4). The trajectories are shown in the x-y coordinates, with a red star at the origin (where the source is). While Figure 2 shows the vehicle successfully learns to reach the stimulus source in 9 out of the 10 trials, the degree to which it approaches the source changes across trials varies. All of the successful simulations stop relatively close to the source. Looking at the figure showing the RBF Network outputs we can see where the movement of the vehicle converges to, by looking at the zero crossings of the function. For instance the RBF network corresponding to the black dashed graph makes the vehicle stop when the stimulus S(x) reaches a stimulus value of 3, and F (s) is negative for larger values of the stimulus, which does not completely fulfil the conditions imposed for the sensori-motor connection function. The RBF network corresponding to the cyan graph, on the other hand, makes the vehicle stop when the stimulus S(x) reaches 3.5, which is much closer to the source. On the other hand, as the figure on the right shows, most of the RBF networks approximate functions with a negative slope, and it is positive for points not experienced during the training. In general, the attractor of the Braitenberg vehicle on the workspace can be obtained by solving S(x) = s 0 , where s 0 is the value at the zero crossing of the RBF network.

In our second set of simulations we used the stimulus sum reward function which penalises perpendicular directions, results of which are shown in Figure 4. For this experiment 8 out of the 10 trials successfully located the stimulus, again with varying degrees of success. However, it should be noted that still most of the RBF networks generate attractors around the maximum, not at the maximum.

In our third set of simulations we used a pose based reward, consisting of the proximity of the vehicle to the source and its heading. Figure 5 shows the results of this simulation. As can be seen from the RBF network, the functions that have been approximated are more accurate, allowing the robot to get much closer to the source before stopping, i.e. before the attractor set of the dynamical system. Using the pose based reward function proved to be advantageous as shown by this set of simulations, with a 100% learning success rate across the 10 trajectories.

In our fourth and final set of simulations we used an altered pose based reward to try to get the vehicle to reach the stimulus source, detailed in Figure 6. As can be seen from the RBF network output, this is the first experiment whereby the RBF networks make the vehicle reach the maximum, i.e. the attractor is at the source. The results of this can be observed with the successful trajectories plotted, where each vehicle easily finds it's way to the stimulus source regardless of where it starts and gets very close to it in the simulated time.

Inverse distance stimulus

Another common example of Braitenberg vehicles is the one implementing phonotaxis through microphones [START_REF] Shaikh | From ear to there: a review of biorobotic models of auditory processing in lizards[END_REF], for instance with a sound source placed at some height h 0 above the ground. It can be seen that, according to the inverse-distance law, the sound intensity will fulfil the conditions for the stimulus function. Furthermore, if the emission pattern of the sound source is isotropic, the stimulus S(x) will be such that S(x) ∝ 1 √ h 2 0 +x 2 +y 2 . This section presents the results obtained for a stimulus source following the inverse distance law, i.e a stimulus of the form S(x) = g0 √ 1+η||x|| 2 . We selected in all the sound-like simulations this last functional form and used g 0 = 4 and η = 0.25. The same experiments were repeated to test whether changing the stimulus affects the performance of the learned controller, i.e. whether some stimulus might be more difficult to follow. As can be seen from the figures below, generally the performance was similar in these sets of simulations compared to the ones above using the inverse-square law stimulus. Like in the phototaxis experiments the RBF networks approximate a function which crosses the zero velocity before reaching the maximum stimulus, meaning the vehicles stop before they reach the stimulus source. The stimulus sum reward function (Figure 7) has the same number of successful trials while the stimulus sum reward function which penalises perpendicular directions (Figure 8) has 10 successful trials, 2 more than the phototaxis counterpart. This proves that, when a sound-like stimulus is used, including a penalty for perpendicular directions helps the Reinforcement Learning algorithm. Interestingly, although the learning functions display a shape different from the previous results (and seem to have a smaller basin of attraction) the connection function F (s) for reward function 3 and 4 generate a point attractor at the stimulus source. Figure 9 shows 10 successful trials using the first pose based reward, based on the proximity of the vehicle to the source and its heading (see Figure 3 for the evolution of this reward), while Figure 10 shows 7 successful trials of the altered pose reward function. Although Figure 10 has less successful trials, the function it approximates is clearly more accurate as the vehicle trajectory is closer to a shortest path trajectory compared to the results that can be seen in Figure 9.

Conclusions and Future Work

This paper presents the first application of Reinforcement Learning to Braitenberg vehicles, a model of insect navigation. Previous robotic works implementing this taxis strategy used hand tuned parameters, or evolutionary strategies [START_REF] Bicho | The dynamic approach to autonomous robotics demonstrated on a low-level vehicle platform[END_REF][START_REF] Mondada | Evolution of neural control structures: some experiments on mobile robots[END_REF][START_REF] Salomon | Evolving and optimising braitenberg vehicles by means of evolution strategies[END_REF] to optimise some cost function dependent on the parameters of a neural controller. We used a reinforcement learning framework to adjust the weights of a radial basis function network to control the vehicle. The presented approach has two main advantages over existing techniques. First, it allows us to define an objective function instead of relying on the perception of the resulting behaviour by the developer. Second, the optimal solution is reached several orders of magnitude faster than the evolutionary approach. On the other hand, the analytic model of Braitenberg vehicles allows to draw conclusions about the stability of the learnt RBF controller.

The results obtained here are limited to computer simulations of the vehicles. In the future we plan to test this framework in real robots performing different types of target reaching behaviours. As an intermediate step, however, we plan to introduce noise in the simulated sensors, i.e. add noise to the stimulus function. This will turn the motion equation of the vehicle into a stochastic differential equation, and theoretical results on deterministic systems might not be applicable. Some early results on the behaviour of vehicle 3a under noise conditions exist [START_REF] Rañó | A drift diffusion model of biological source seeking for mobile robots[END_REF] that can help interpreting the results of (or simplifying) the learning problem.

Fig. 1 .

 1 Fig. 1. Braitenberg vehicle 3a.

I0 h 2 0

 2 +||x||2 , where I 0 is the intensity of the light at the source. This function can be rewritten as S(x) = g0 1+η||x|| 2 where η = 1 in all the light-like simulations this last functional form and used g 0 = 4 and η = 0.25.

Fig. 2 .

 2 Fig. 2. Results from phototaxis simulations using reward function 1.

Fig. 3 .

 3 Fig.3. On the left, an example from the phototaxis simulations of the evolution of a stimulus based reward (using reward function 1), and on the right, an example from the phonotaxis simulations of the evolution of a pose based reward (using reward function 3).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Results from phototaxis simulations using reward function 2.

Fig. 6 .

 6 Fig.6. Results from phototaxis simulations using reward function 4.

Fig. 7 .

 7 Fig. 7. Results from phonotaxis simulations using reward function 1.

Fig. 8 .

 8 Fig. 8. Results from phonotaxis simulations using reward function 2.

Fig. 9 .

 9 Fig. 9. Results from phonotaxis simulations using reward function 3.

Fig. 10 .

 10 Fig. 10. Results from phonotaxis simulations using reward function 4.

Acknowledgements

This work was partially supported by the Royal Society International Exchange Scheme under grant IE151293.