A review of laboratory and numerical modelling in volcanology

Modelling has been used in the study of volcanic systems for more than 100 years, building upon the approach first applied by Sir James Hall in 1815. Informed by observations of volcanological phenomena in nature, including eyewitness accounts of eruptions, geophysical or geodetic monitoring of active volcanoes, and geological analysis of ancient deposits, laboratory and numerical models have been used to describe and quantify volcanic and magmatic processes that span orders of magnitudes of time and space. We review the use of laboratory and numerical modelling in volcanological research, focussing on sub-surface and eruptive processes including the accretion and evolution of magma chambers, the propagation of sheet intrusions, the development of volcanic flows (lava flows, pyroclastic density currents, and lahars), volcanic plume formation, and ash dispersal. When first introduced into volcanology, laboratory experiments and numerical simulations marked a transition in approach from broadly qualitative to increasingly quantitative research. These methods are now widely used in volcanology to describe the physical and chemical behaviours that govern volcanic and magmatic systems. Creating simplified models of highly dynamical systems enables volcanologists to simulate and potentially predict the nature and impact of future eruptions. These tools have provided significant insights into many aspects of the volcanic plumbing system and eruptive processes. The largest scientific advances in volcanology have come from a multidisciplinary approach, applying developments in diverse fields such as engineering and computer science to study magmatic and volcanic phenomena. A global effort in the integration of laboratory and numerical volcano modelling is now required to tackle key problems in volcanology and points towards the importance of benchmarking exercises and the need for protocols to be developed so that models are routinely tested against “real world” data.


Introduction
Many people have a profound fascination with volcanoes: from the beautiful landscapes they produce and art they inspire (Sigurdsson, 2015b) to the impact their eruptions have on individuals, societies, and civilisations (e.g.Sheets, 2015).However, this fascination is often met with fear as the destruction volcanoes can cause may have far-reaching effects.Volcanic activity is often unpredictable and occurs in an environment that is highly changeable and forbidding; however, there is a compelling need to improve our understanding of these complex systems.Approximately 800 million people around the world live close enough to a volcano to be directly affected by an eruption (Loughlin et al., 2015), and many more are at risk of social or economic impact as the consequences of volcanism extend from regional to potentially global reaches (e.g.Svensen et al., 2004).The huge range of style and intensity of volcanic activity means that societies living nearby do so with high risk; there are however many benefits too: volcanoes produce habitable environments on a local scale through production of fertile soil, and magmatic activity is associated with economic deposits such as copper porphyry.On a planetary scale, the gases volcanoes emit lead to the creation of our oceans and the atmosphere.Life on Earth and the physical processes that govern volcanic activity are thus intimately connected.
The effects of volcanic eruptions can be felt long after surface activity has ceased, with the potential for both physical environments and societies to be impacted many decades after the event (e.g. the occurrence of lahars decades after the 1991 Mt Pinatubo eruption).The challenges of working in volcanic terrains and gathering useful data mean that laboratory and numerical models have gained significant importance in studying the dynamics of volcano growth and erup-Published by Copernicus Publications on behalf of the European Geosciences Union.
tion.The occurrence of volcanic and magmatic activity is challenging and potentially impossible to forecast, and the factors that influence processes such as magma storage, ascent, eruption, and deposition have several, often poorly constrained, variables.Technological limitations place boundaries on what we can record in nature; for example, direct observations are often problematic as the processes of interest are frequently hidden from view as they occur beneath the Earth's surface or within a pyroclastic flow or plume.The often remote and difficult-to-access location of volcanoes poses further logistical problems.
Volcanic and magmatic processes occupy a vast range of scales from sub-millimetre to kilometres in size, and the timescales over which they take place range over fractions of a second up to millennia.The benefits of laboratory and numerical modelling mean that it is possible to study these processes in a controlled way, with the opportunity to repeat the phenomenon as required.Laboratory experiments and numerical simulations can link well-constrained starting conditions with measurable outcomes, and with careful scaling these findings can be extrapolated to better understand the natural processes.
In this invited review, we summarise laboratory and numerical modelling in the broad field of volcanology by providing some historical context and overview of models of magma and lava rheology, magma chambers, magma intrusions, lava lakes and lava domes, lava flows, pyroclastic density currents, lahars, volcanic plumes, and ash dispersal (see Fig. 1).Our paper is aimed at graduate students and those with a general interest in modelling volcanoes.Our intention is not to consider all parts of the magmatic and volcanic system or review each area in great depth.Instead we provide an account of the foundations of the subject, highlight some key papers, and describe some of the cutting-edge techniques that are being utilised to model volcanic and magmatic processes numerically and in the laboratory.We hope that this contribution inspires future reviews on topics we have touched on or been unable to include, and particularly with a focus on developments in both numerical and laboratory modelling.We conclude by highlighting common challenges in the future of volcanology in our efforts to model these intriguing, captivating, and potentially devastating phenomena.
2 Historical context of modelling in volcanology 2.1 Volcanology as a science Our fear and fascination with volcanoes is evident in some of the earliest historical accounts of volcanic eruptions, which show that we have long tried to understand how volcanoes form and what causes them to erupt.Greek natural philosophers from the fifth and fourth century BC, such as , , proposed that volcanic eruptions Figure 1.Schematic illustration of the magmatic and volcanic phenomena that have been modelled in the laboratory and in numerical simulations.In this review we focus on (I) magma chambers (see Sect. 5), magma sheet intrusions such as (II) sills and (III) dykes (see Sect. 6), lava lakes (Sect.7), lava domes (see Sect. 8), (IV) lava flows (Sect.9), (V) pyroclastic density currents and lahars (see Sect. 10), (VI) volcanic plumes (see Sect. 11), and (VII) volcanic ash dispersal (see Sect. 12).
were caused by "great winds inside the Earth", an idea that was supported by Aristotle (384-322 BC) (see Sigurdsson, 2015a, for details).The study of volcanoes is founded on eye-witness accounts and field observations of activity, beginning in Italy with the descriptions of Pliny the Younger of the eruption of Vesuvius in 79 AD.Centuries later, William Hamilton (a pioneer in volcanology) gained recognition for his descriptions of the eruption of Vesuvius in 1767 (Hamilton and Cadell, 1774).The interpretation of the origin of intrusive and extrusive igneous materials in the rock record came shortly after with the seminal work of James Hutton, the so-called "Father of Geology", whose field observations in 1785 in the Cairngorm mountains, Scottish Highlands, demonstrated that granite was formed by the solidification of an initially liquid body that had intruded a pre-existing host rock (Hutton, 1788).

Experiments in volcanology and a quantitative approach
Experimentation has long been an important tool in geological investigations.The first analogue experiments to study geological processes were published just over 200 years ago by James Hall in 1815.Hall's experiments, carried out "with such materials as were at hand", used several pieces of cloth, linen and wool, each placed upon one another horizontally to represent a sequence of rock strata (see Fig. 2).A rig then applied horizontal shortening to create a folded structure that Hall deemed reminiscent of convoluted rock layers at Fast Castle, Cockburnspath, Scotland.Hall's research was read at a meeting of the Royal Society of Edinburgh, and has since inspired a new field of study using analogue materials to explore geological processes in the laboratory.Several decades later, the eminent experimentalist Auguste Daubrée used a device to study the formation of cavities in cylindrical rock samples due to the explosive release of gas, using his results to explain the development of volcanic diatremes related to diamondiferous kimberlite deposits (Daubrée, 1891).
He used dynamite or guncotton as analogues for the rapid release of volatiles from ascending magma or explosive energy due to interaction with water, and so Daubrée's experiments laid the foundations for the use of analogue laboratory modelling in volcanology.
It was not until the 1960s that quantitative science started to emerge more prominently in the volcanological literature, transitioning away from what had been largely qualitative and descriptive work.The use of numerical and analogue modelling in volcanology has come from the application of methodologies developed for alternative, sometimes quite disparate, purposes.From the 1950s the fields of volcanology, fluid dynamics, and engineering came together (Fig. 3).Hubbert and Willis (1957) used gelatine solids injected with plaster of Paris to study hydraulic fractures formed within pressurised boreholes; however, the experiments also apply to modelling magma-filled fractures such as dykes and sills (Fig. 3a).Fiske and Jackson (1972) subsequently used freestanding gelatine models to study magma propagation in a volcanic rift such as Hawaii (Fig. 3b).The research of Morton et al. (1956) modelling industrial plumes rising from chimney stacks is the basis of many of the numerical models of volcanic plumes implemented today (Fig. 3c).
A pioneer in quantitative approaches in volcanology was George Walker (1926Walker ( -2005)), who is considered by many in the field to be "the father of modern quantitative volcanology" due to his demonstration of how to integrate disciplines across the sciences, using field-based measurements to develop and test conceptual models.The breadth of Walker's interests and expertise across volcanology were impressive, and his guidance and influence on many parts of the volcanic system from plumbing systems, lava flows, tephra fall, and flows is evident in the literature (see Sparks, 2009, for a review).Walker's quantitative approach, for example mapping  Hall (1815).A series of images depict a set of experiments which have since inspired the use of analogue materials to provide a physical explanation for geological observations made in the field.Hall produced an "ideal" coastal section (a) to demonstrate the continuous nature of the folded rock layers observed in the field.Sketches depict a set of experiments that were performed (b) by compressing clay layers to produce convolutions (c) that are reminiscent of the structures observed in the field.
the regional distribution of zeolites in Icelandic basalt lavas (Walker, 1960) and being one of the first to publish calculations of lava flow viscosity based on field measurements of flow thickness, velocity, and angle of slope (Walker, 1967), led to the rejuvenation of laboratory experiments in volcanology (and eventually numerical simulations) that could test hypotheses based on field data.George Walker and Lionel Wilson published laboratory experiments that studied the physics of pyroclastic fallout from large and highly explosive volcanic eruptions, focussing on quantifying the rate of fall of ejected debris.They timed the fall of carefully characterized natural samples of tephra (pumice, lithic, and crystals) and compared these with analytical computed terminal velocities of cylinders and rough and smooth spheres to aid the analysis and interpretation of field deposits (e.g.Walker et al., 1971).Steve Sparks and Lionel Wilson developed early analytical models to explore the controls on volcanic column height (Wilson et al., 1978) and the role of vent geometry in the collapse of eruptive columns to form ignimbrite deposits (Sparks and Wilson, 1976).These approaches, of comparing analytical and empirical models with field data, underpin  (Hubbert and Willis, 1957), (b) water injected into a freestanding triangular prism of gelatine to simulate magma intrusion in a volcanic rift (Fiske and Jackson, 1972), and (c) injection of lowdensity fluid into a uniform ambient fluid (Morton et al., 1956).much of the numerical and laboratory modelling work conducted in volcanological research today.

Parameterisation and scaling of models in volcanology
Modelling is a process that searches for the universal laws that control a natural phenomenon.The aim is not to reproduce a specific natural example, but to understand the role of a limited number of parameters on the processes being studied.Analytical models, numerical simulations, and laboratory experiments have helped to explain processes that are too large or too complex to be understood in nature.They provide physical intuition and insight that span a large parameter space by establishing trends, influences, and interdependencies of key parameters.We define a numerical model as a set of algorithms and equations that are used to capture the physical or chemical behaviour of the system being modelled.We define an analogue model as a simplified representation of physical processes that is scaled down so it can be studied in the laboratory.The term "analogue modelling" describes the process of laboratory experimentation using prototype materials.This process-oriented approach aims to identify underlying physical processes that control the studied phenomenon, and as such is sometimes called "laboratory modelling" or "analogue laboratory experimentation".However, this definition can be expanded to include recent developments that have enabled dramatic large-scale experiments using natural materials (e.g.Lev et al., 2012).Recent advancements in computational power, analytical techniques, and experiment imaging have revolutionised these methods, ensuring their continued use in studying volcanic and magmatic phenomena.Volcanic and magmatic processes are controlled by a range of physical processes and regimes, and these vary depending on whether subsurface or eruptive processes are being con-sidered.For example, subsurface processes such as magma intrusion and conduit processes are largely controlled by the rheology of magma and deformation of the host rock, while eruptive processes are dependent on the characteristics of the eruptive mixture and its relation to the ambient fluid.Modelling volcanic and magmatic mixtures is non-trivial due to their multiphase nature; magmas comprise melt, crystals, and gas, while eruptive plumes and flows commonly comprise many different particles sizes, in addition to a gas phase.Fundamental to the success of models in volcanology is their scaling and parameterisation, with the model outputs strongly dependent on the quality of the model inputs.However, there are different considerations that need to be made when developing a numerical model or laboratory experiment.

Scaling
Scaling is important in both laboratory modelling and numerical modelling.Ampferer (1906) was the first to recognise that when geometric lengths are reduced so that a phenomena can be studied in the laboratory, then all other "constants" should also be reduced.The fundamental principles of scaling are that physical phenomenon can be described by equations, and that these equations must be balanced dimensionally (their units are homogeneous).The fundamental units are length [L] in metres, mass [M] in kilograms, and time [T ] in seconds (where square brackets mean "with units of").Additional fundamental units may include temperature [K] in degrees Kelvin.These fundamental units combine to produce many derived units, such as velocity [LT −1 ], acceleration [LT −2 ], pressure (or stress) [ML −1 T −2 ] (or Pa), viscosity [ML −1 T −1 ] (or Pa•s), and energy [ML 2 T −2 ] (or J ).
Scaling procedures are different depending on the part of the volcanic system that is being studied, but there are two main approaches that are used: (1) scale factors and (2) dimensional analysis (also known as Pi theorem).A brief description of these methods is provided below, but the reader is directed to Merle (2015) and Galland et al. (2018) for more detailed accounts and descriptions of their application in volcanology (both for numerical and analogue modelling).

Scale factors
Scale factors have been applied to a range of volcanic and magmatic phenomena, considering the ratio of characteristic parameters in nature and the model.The model ratios describe the conditions for geometric, kinematic, and dynamic similarity, and can be defined by a combination of a length ratio h * , time ratio t * , and mass ratio m * : where lower case "m" and "n" indicate the model and nature, respectively.Geometric similarity means that the ratio of distances in the model and nature is a constant, and kinematic similarity means that this geometric similarity is maintained over time.The time ratio then remains constant for the whole experiment (e.g. the duration of an experiment scales with the duration of the process in nature).Dynamic similarity includes a wide range of force ratios F * , such as gravity or viscosity.Directions in the model and nature need to be the same (for example when scaling velocities or mechanical forces).

Dimensional analysis
Originally described by Buckingham (1914), dimensional analysis (also known as Pi theorem) assesses similarity of a model and its natural counterpart by comparing a set of dimensionless numbers.It can be used to identify physical regimes in volcanological problems.Pi numbers ( ) are defined as a set of independent dimensionless numbers.They are constructed from a set of k independent dimensional numbers that are based on i variables (or parameters) and j fundamental units (k = i − j ).A commonly used Pi number is the Reynolds number: Here ρ is fluid density, u is the flow velocity, L is the characteristic length, and µ is the viscosity.It describes the ratio of inertial and viscous forces, defining laminar flow at low Re, turbulent flow at high Re, and a transitional regime between.The Rayleigh number Ra is also a frequently used Pi number in volcanology: where g is gravity, α is the coefficient of thermal expansion, κ is the thermal diffusivity, and T is the temperature difference.It describes the ratio of buoyancy and viscous forces, and quantifies the transition from heat loss by conduction to heat loss by convection.For a full description of Pi numbers and how they are constructed, the reader is referred to the original text by Buckingham (1914) or recent explanation by Merle (2015).Similarity is ensured when the calculated Pi numbers are identical or very similar in the model and in nature, even if the values in the governing dimensional parameters are very different.However, if the individual Pi numbers do not agree, then the model and nature are not similar; but it can still be argued that the processes are physically equivalent if when plotted against one other they lie on the same scaling law (e.g.Galland et al., 2018).The most challenging aspects of dimensional analysis and using Pi theorem are selecting the variables that are important in the phenomenon being studied and determining which are the most useful dimensional numbers they produce.

Numerical modelling approaches
Numerical modelling involves the construction of a number of mathematical, physical, and/or chemical equations and simplifications to represent a particular phenomenon (Barenblatt, 2003).These equations are then solved according to a set of initial and boundary conditions, and they are modified or developed based on the model assumptions.
Deterministic models aim to reproduce the underlying physical behaviour of a process, where the model output is controlled by the model parameters and the initial conditions.Application of deterministic models is difficult due to the simplifications that need to be made due to the complexity and limited knowledge regarding volcanic systems, for example the use of coefficients to describe entrainment of ambient air into a rising volcanic plume.Stochastic models are applied to capture the probability of various outputs and incorporate an inherent randomness such that the same set of parameter values and initial conditions will produce a range of model outputs.In volcanology, most numerical models are deterministic and the application of stochastic approaches is generally applied to hazard assessments.However, application of statistical modelling approaches such as Monte Carlo or Latin hypercube simulations enable deterministic models to be used in model uncertainty quantification or sensitivity analysis.These methods investigate the variance in model results when inputs are changed, identifying parameter dependencies and those parameters that have a first-order control on model outputs.
Numerical models with a range of complexity exist and have different uses.Steady-state models provide a first approximation of the physical behaviour and a simplified view of the processes involved, requiring a number of assumptions to be made.In comparison, transient models consider system changes over time, are more complex, and require the closure of more equations, which makes them more computationally expensive.However, they are able to reproduce the unsteady behaviour observed in a range of volcanic phenomena.
Regardless of complexity, all numerical models require description of boundary conditions which take into account interaction with the surrounding environment.Examples of boundary conditions relevant to volcanic systems include atmospheric conditions when modelling eruption plumes, as humidity and wind strength exert a strong control over the height a volcanic plume may reach in the atmosphere, and stiffness of a host rock and magma viscosity when modelling magma intrusions, as these will affect the propagation behaviour of a dyke or sill.Environmental conditions can be approximated in simplified numerical models, for example when modelling magma chamber growth the lithosphere may be modelled as a mechanically homogeneous material (e.g.Galgana et al., 2013) or by including mechanical heterogeneity and a stiffness contrast between crust and mantle (e.g.Le Corvec et al., 2015).Interactions between a phenomenon and its surroundings are commonly accounted for in numerical models by simplified coefficients; for example, the entrainment of ambient air into a rising plume is described by two entrainment parameters: one considers radial entrainment due to turbulent eddies at the plume edge, and the other accounts for entrainment due to the effects of wind on the plume.These coefficients have been parameterised using a combination of observations of the phenomena in nature, analogous phenomena such as chimney stacks (Morton et al., 1956), and laboratory experiments (Fig. 4).While it is possible to account for the effect of the ambient conditions on the modelled phenomenon, it is considerably more difficult to account for feedback between the modelled phenomenon and its environment as this requires significantly more computational power.

Laboratory modelling methods
The choice of analogue materials for laboratory modelling will depend on the parameters that are being investigated and the conceptual model that is being tested; many simplifications are required in order to track the impact of variables on experiment outcomes, e.g.how changes in density contrast between fluid and surroundings effect geometry, velocity, or pressure.Extracting and measuring parameters and variables in laboratory experiments is crucial for understanding the modelled phenomenon, and for checking model parameterisation is consistent with appropriate scaling laws.Prior to this, careful material characterisation is required at experimental conditions.

Analogue materials and experiment imaging techniques
There have been several recent developments in imaging and measuring experiment materials and variables in the laboratory, drawing on technologies developed for industrial purposes to study volcanic processes.Imaging and measurements need to focus on detailed external and internal monitoring of a laboratory experiment.
To assist with analogue material selection, temperaturedependant Newtonian and non-Newtonian fluids and gels are characterised in a viscometer or rheometer at a controlled temperature and using a range of measurement geometries; these include rotary, concentric cylinder, falling ball, tube A flow diagram to represent the optimal approach to use analogue and numerical modelling in volcanology.Observations and measurements from natural volcanic and magmatic phenomena provide the parameters and initial conditions for analogue and numerical models; these models are then tested against nature, with analogue models also aiding the development of numerical models.and parallel plate methods.Material properties can often be tailored so that it meets the required physical behaviour.
Photogrammetry enables measurements to be made from photographs with great precision and can be used to construct three-dimensional representations of real-world objects.Structure-from-motion algorithms process experiment images from synchronised cameras and, for example, create a time series of digital elevation models representing the changes in the surface elevation of an experiment.
X-ray micro-tomography can produce high-resolution three-dimensional reconstructions of a static model topography with cross sections showing internal structures.
Techniques such as digital image correlation (DIC) can map strain changes and particle image velocimetry (PIV) or optical flow techniques can be used to map fluid flow; lasers can also be used to illuminate a thin vertical sheet within a gelatine experiment or particle suspension experiment, with images recorded at time-defined intervals.
The ability to image and integrate measurements of the surface and subsurface development of, for example, magma intrusion in the laboratory greatly strengthens the ability of analogue experimentation to help inform numerical modelling that is used to interpret volcano-deformation data in nature.

Magma rheology
The methods described in Sect. 3 have been used to study magmatic and volcanic processes that extend from the depths of the lithosphere into the stratosphere.In the following sections, some of the numerical and scaled laboratory models of magma chambers, sheet intrusions, lava lakes, lava domes, lava flows, pyroclastic density currents, lahars, vol-canic plumes, and ash dispersal are described.Magma is one of the principal components of all volcanic systems, and how it behaves as it flows and fragments is key to tackling perhaps all processes in volcanology.Therefore, our review begins with a detailed account of modelling magma rheology.
Magma is a multiphase fluid, comprising a melt phase with variable proportions of bubbles and crystals.The viscosity (η) of pure melt is considered Newtonian (Fig. 5), with a linear relationship between shear stress (σ ) and strain rate ( γ ) (Lejeune and Richet, 1995;Ishibashi, 2009): Due to its multiphase nature, magma is considered non-Newtonian.Several types of non-Newtonian rheology have been applied to model the behaviour of magmas and lavas based on field observations (Fig. 5).A Bingham fluid has to overcome a yield stress before it can begin to flow (Hulme, 1974): where σ 0 is the initial shear stress required to cause the onset of flow when γ = 0. Once the yield stress has been overcome, the fluid has a constant viscosity.More recently, the Herschel-Bulkley model (Herschel and Bulkley, 1926) has been applied to the behaviour of magmas due to its versatility in allowing for the modelling of a spectrum of magma behaviours (Llewellin et al., 2002b;Mueller et al., 2011): where σ 0 is yield stress when there is no flow, K is the consistency (η when γ = 1), and n is the degree of non-Newtonian behaviour (where n = 1 is Newtonian, n<1 is shear thinning, and n>1 is shear thickening).Viscosity is a key parameter that is considered in several of the important dimensionless numbers used in scaling experiments in volcanology, such as the Reynolds number, Rayleigh number, and Peclet number.Depending on the application and level of complexity, a variety of analogue materials have been used to model magma (see Table 1 for a summary).An important online resource administered by Alison Rust (University of Bristol) catalogues a range of analogue materials that have been used in volcanology, describing their material properties, applications, and limitations (https://sites.google.com/site/volcanologyanalogues/home).Many models use a meltonly magma analogue for simplicity, or in more complex models two-phase suspensions (bubbles in liquid, or crystals in liquid) and rarely three-phase (bubbles and crystals in a liquid).As such, the spectrum of rheology that has been considered in magma analogue models is broad and includes the use of Newtonian fluids, Bingham fluids, or Herschel-Bulkley fluids.4.1 Two-phase suspensions: particle suspensions Particulate suspensions are ubiquitous across the volcanic system; from crystals in magma, where variations in crystal content mostly originate from changes in temperature, to ash particles within an eruptive plume, where the particle size and concentration within a volcanic plume or pyroclastic density current is related to the type of eruption.For the purposes of numerical and laboratory modelling, particle distributions are simplified by a single well-defined particle size or a small number of particles sizes used to replicate natural systems (e.g.Fig. 6a-c).However, both numerical and laboratory studies have shown that particulate concentration has a first-order control on eruptive behaviour.The Einstein-Roscoe equation has been used in magma modelling to describe how the increase in crystal content affects the bulk viscosity, where an increase in the particle volume fraction (φ p ) causes an increase in bulk magma viscosity.Low φ p magmas behave as Newtonian fluids whereas high φ p magmas (φ p = >30−40 %) behave as non-Newtonian fluids (e.g.Lejeune and Richet, 1995;Mueller et al., 2011).Bimodal crystal size populations act to reduce melt viscosity, with higher viscosities observed in unimodal crystal size suspensions compared with lower viscosities in bimodal crystal size suspensions (Castruccio et al., 2010).In volcanic plumes, high particulate concentrations, leading to high plume densities relative to the atmosphere result in column collapse (e.g.Carey et al., 1988).Both numerical and analogue studies have shown that particle concentration has a first-order control on initiation of coignimbrite plumes (Andrews and Manga, 2011;Engwell et al., 2016).The settling behaviour of particles has a significant control on the behaviour of propagating currents, particularly ash plumes and pyroclastic density.For example, within ash plumes, larger particles are deposited closer to the source, while smaller particles are re-  Truby et al., 2015).(e) Bubble injection experiments using a small-gap parallel-plate geometry to study the development of permeable pathways in a particle-rich suspension.Particle image velocimetry has been used to measure particle speed in three experiments with different crystal fraction (Oppenheimer et al., 2015).tained over much greater distances, controlled by the particle terminal velocity where g is gravitational acceleration, ρ p and ρ g are the densities of the particles and the surrounding gas, respectively, d is the particle diameter, and C d is the drag coefficient determined from analogue experiments (Kunii and Levenspiel, 1991).The terminal velocity can be obtained from an experimental correlation of the dimensionless groups C d vs. the particle Reynolds number, Re, defined as where µ (Pa•s) is the dynamic viscosity of the gas.For spherical particles, a number of analytical equations can be used to determine C d depending on the particle Richardson number (Bonadonna and Phillips, 2003).In particulate suspensions, the concentration of particles also acts to control settling behaviour, with high concentrations leading to the development of Rayleigh-Taylor, or gravitational instabilities, as those regions become denser than the ambient air, as shown by a number of analogue laboratory studies (Carey et al., 1988;Carazzo and Jellinek, 2012;Manzella et al., 2015).

Two-phase suspensions: bubble suspensions
The effect that bubbles have on magma viscosity depends on bubble shape, size, and ability to deform under stress.In steady flow regimes, where stress and shear are constant, the bubbles reach an equilibrium deformation that can be measured by the dimensionless capillary number Ca, which describes the relative effects of viscous forces and surface tension acting between a liquid and a gas (Manga and Loewenberg, 2001;Llewellin et al., 2002b): where η 0 is the fluid viscosity without bubbles, r is the undeformed bubble radius, γ is strain rate, and is interface surface tension between the liquid and gas.Small capillary numbers are dominated by surface tension, meaning that bubbles reach their equilibrium deformation soon after there is a change in shear rate (Llewellin et al., 2002a, b).This produces spherical bubbles that act to increase the viscosity of the suspension by creating an obstacle to flow.Large capillary numbers give rise to easily deformable and often elongate bubbles, acting as sites where shear localisation can occur due to a reduction in friction, and therefore the presence of large bubbles will have a reduced impact on bulk viscosity (Manga et al., 1998;Mader et al., 2013).
In unsteady flow regimes, when there is a variable strain rate, the forces causing deformation and restoration of the bubble shape are not in equilibrium (Llewellin et al., 2002b).As such, the Ca number (Eq.11) does not adequately describe the behaviour of the bubble, and so a dynamic capillary number Cd is defined: where γ is the rate of change of the imposed deforming force and is the relaxation time, which is the time taken for the bubble to return to spherical under the action of surface tension (Mader et al., 2013):

Three-phase suspensions
Three-phase suspensions are well suited to explaining the behaviour of magmas and bring us closer to understanding volcanic systems, but they also present several challenges associated with the additional complexity modelled.A threephase suspension can be modelled assuming a bubble suspension base fluid with the particles suspended within (Truby et al., 2015; see Fig. 6d): where η * is relative viscosity, η b is bubble suspension viscosity, ϕ p is particle volume fraction and ϕ m is the maximum packing fraction.This simplifies the calculation of the three-phase rheology and assumes a low Ca (see Truby et al., 2015); however, if this is not appropriate then η b may need to be substituted by either the high Ca or Cd viscosity equation.This new model can account for a crystal-bearing magma that has no bubble content at depth but vesiculates during ascent.
It therefore marks a significant advancement in our understanding of magma behaviour through time and space and will be an important tool in future models to better constrain the impact of three-phase magma rheology on volcanic eruptions.Laboratory experiments of degassing crystal-rich magmas have shown using analogue materials that gas escape and the development of permeable pathways in particle-rich suspensions can be fracture-like or due to bubble formation, and that the migration pathways are controlled by particle fraction and the degree of particle packing (see Fig. 6e; Oppenheimer et al., 2015).

Magma chambers
The largest accumulations of magma on Earth are thought to occur in deep reservoirs called magma chambers that have the potential to feed volcanic eruptions or stall to form a pluton (e.g.Lundstrom and Glazner, 2016a;Purtirka and Cooper, 2017).Magma chambers are perhaps the most obscure parts of the volcanic and magmatic system.Scaled models have considered heat and flow conditions, focussing on particular regimes by assessing the impact of dimensionless numbers such as the Reynolds number, Rayleigh number, and Peclet number on magma chamber processes.Plutons are large accumulations of coarsely grained igneous rock and they express a broad range of compositions, generally falling between granite and gabbro, with physical properties (such as viscosity) that can span several orders of magnitude and vary in both space and time.The relationship between large, ancient magma bodies, such as laccoliths and plutons, and magma chambers that feed volcanic eruptions is enigmatic and currently under debate (e.g.Lundstrom and Glazner, 2016b).A growing body of literature challenges the traditional "big tank" conceptual model of magma chambers as dynamic, large, and long-lived accumulation of magma that slowly cools, crystallises, and differentiates (Glazner et al., 2004).The relevance of experiments that study pluton emplacement as a single body of deformable, viscous fluid in the crust therefore needs to be reassessed.Instead it is proposed that large igneous bodies and magma chambers are incrementally emplaced from the accumulation of sill-like bodies (horizontal planar magma-filled sheets), and that they are discrete and ephemeral regions of melt and mush (see Annen et al., 2015, for a review).

Building magma chambers
Numerical models have explored the wide parameter space that permits the existence of magma chambers and their ability to feed volcanic eruptions, aiding the interpretation of the origin of plutons in the rock record but also the deformation currently manifested at large volcanic systems.The application of magma chamber numerical models ranges from constraining how the Earth's crust is built to explaining how magma chambers can be generated to feed the largest eruptions on Earth (e.g.Lundstrom and Galzner, 2016a).Thermomechanical numerical models study the physical processes that build magma chambers, for example by the focussing of dykes towards a magma chamber.The conceptual model being explored is that crustal magma chambers are built incrementally by repeated injections of dykes and sills.Karlstrom et al. ( 2010) modelled a magma chamber fed by dykes from below that intruded into the "capture radius" of the chamber to grow the magma chamber.They emphasised the two-way coupling between magma chamber stress and thermal evolution by modelling the magma chamber as a circular fluid-filled cavity surrounded by a viscoelastic shell in an infinite elastic material.Karakas and Dufek (2015) developed a thermo-mechanical model that studies the relative roles of magma flux and tectonics in building magma reservoirs in the crust, finding that magmatism and tectonics act together to modify the thermal, mechanical, and chemical structure of the crust.They found that magma flux has a primary control on melt generation, but that an extensional environment helped to amplify this.Karakas et al. (2017) used a finite volume scheme model to simulate a 2-D crustal section 80 km wide and 60 km deep with 40 m resolution.They constrained conditions for lower and upper crustal reservoirs and how the presence of a lower crustal magma reservoir modifies the thermal structure of the crust and acts to reduce the magma flux required to sustain upper crustal magma chambers.Complementary to thermo-mechanical numerical models are thermodynamic numerical models, which provide equilibrium constraints on the petrological and geochemical evolution of the system.The "magma chamber simulator" of Bohrson et al. (2014) is a sophisticated thermodynamic numerical model that calculates phase equilibria and geochemical evolution of resident magma, recharge, and wall rock in a self-consistent way.Simultaneous magma recharge, assimilation, and crystallisation are considered in a system comprising melt + solid + / − fluid, and the model outputs can help constrain the interpretation of natural geochemical datasets.A series of papers by Annen et al. (e.g. Annen et al., 2006a, b;Annen, 2009) account for the generation of intermediate and silicic magmas in deep crustal hot zones (Deep Hot Zones), based on a heat-transfer numerical model that simulates magma injection in the crust or crust-mantle boundary and calculates the conditions required for the accumulation of melt to build a reservoir of eruptible magma.Their equilibrium thermodynamic numerical model is based on a heat balance between injected magma and surrounding rock, and parameters such as density, specific heat capacity, temperature, time, melt fraction, latent heat of fusion, thermal conductivity, and depth are included.These numerical models have identified magma flux as a key parameter that enables a magma chamber to exist, and they calculated that a 10 km thick pluton requires a magma flux that exceeds 10 −2 km 3 yr −1 to permit the development of a magma chamber with a volume of eruptible magma sufficient to feed the largest silicic explosive eruptions.

Magma chamber dynamics and recharge
Several papers were published in the 1980s exploring the cooling of a large predominantly liquid magma reservoir using laboratory models, for example studying the so-called "double diffusive convection" model and its application to magma chamber evolution (Huppert and Turner, 1981b).The double diffusive convection model accounts for convection of a fluid with composition and temperature gradients acting in opposing directions.It was originally derived from oceanographic applications, yet has proved helpful to explain geological features such as large-scale cyclic crystal layering within large igneous bodies and supports the proposed conceptual model that magma chambers are compositionally zoned (Huppert and Sparks, 1980).Subsequent laboratory models have investigated magma mixing and magma mingling, for example for the case of a rhyolite magma chamber injected by basalt from below (Huppert et al., 1983) or mafic magma chambers replenished by felsic injections (Weinberg and Leitch, 1998).The impact of volatile exsolution and bubble formation on magma mixing was explored in the laboratory by Huppert et al. (1982) and Turner et al. (1983) by introducing reactive HNO 3 into liquid layers of K 2 CO 3 (upper layer) and KNO 3 (lower layer) to cause the release of gas.Phillips and Woods (2001) then studied the accumulation of bubbles and movement of bubbles within a magma chamber, using a salt solution as the magma analogue and an electrolysis cell with gauze to produce the bubbles.All these studies demonstrated that a recharge event of bubble-rich and lowdensity magma, such as basalt, into a magma chamber may generate a turbulent bubble plume within the chamber that can be described by plume theory.Such bubble plumes could impact magma mixing within the chamber and the stability of the magma chamber; they could also have the potential to trigger an eruption and could affect the style of eruptive activity.Magma contamination from roof and wall melting has also been studied experimentally (Leitch, 2004).
Over the past decade there has been strong advancement in numerical models of magma chambers that are dominated by melt, due to the incorporation of micro-scale chemical and physical processes.Bergantz (2000) studied a meltrich system using scaling relations to study magma mixing following a magma recharge event.Depending on the Reynolds number (Re), different mixing regimes were iden-tified.At high Re>100 internal contacts between resident and recharge magma were found to be unstable and collapse, causing chamber-wide magma mixing and mingling.Under these conditions, mushy regions within the chamber could be scoured away and previously crystallised material widely distributed.At intermediate Re values (10−100), mixing was incomplete and "islands" of unmixed material remained intact.At low Re values (<1) internal slumping and folding of the magma was expected and would be most likely to preserve a mineral fabric in the rock record.This model did not consider heat transfer, but focussed on the physical effects of heat through changes in viscosity and density.The effects of heat and its impact on the inception of convection during cooling has since been studied (Huber et al., 2008), constraining the timescales for convective overturn and for complete magma mixing and homogenisation during crystallisation (Huber et al., 2009).Ruprecht et al. (2008) considered the impact of gas exsolution on mixing in magma chambers, as gases exsolved from a recharging magma could cause density inversion and overturn of the chamber.They found that such processes could be recorded in the crystal population, with crystallisation recording the initial and final environments but dissolution occurring in the overturn events.

Magma chambers feeding eruptions
Fracture and failure of a magma chamber to feed a volcanic eruption has also been studied in laboratory experiments, studying the nucleation of magma-filled fractures (dykes).McLeod and Tait (1999) used gelatine models to study the pressurisation and failure of liquid-filled cavities, creating a crustal magma chamber by inflating a balloon within the liquid gelatine and removing it when the gel had solidified.Fluid was then injected into the cavity using a head pressure.Under increasing stress and strain, the gelatine undergoes an initially elastic deformation and then brittle failure.They found that dyke nucleation occurred from a pre-existing flaw in the analogue chamber wall, and that the viscosity of the fluid influenced the tendency for dykes to propagate.Canon-Tapia and Merle (2006) carried out similar experiments but injected silicone into a cavity in the gelatine, finding that a sustained overpressure in the magma chamber was required to create dykes that could potentially reach the surface.Koyaguchi and Takada (1994) also used gelatine models and glycerine to explore how the evacuation of a low-viscosity fluid may lubricate the path of a more viscous fluid into a pre-existing fracture.
Earthquakes are potentially an important external influence on magma chamber stability that may trigger dyking that leads to a volcanic eruption.Namiki et al. ( 2016) considered two different scenarios of foam stability over a liquid layer of diluted glucose syrup in a partially filled tank (open vent) or fully filled tank with density-stratified fluids (sealed magma reservoir).The use of a shaking table enabled the authors to identify the conditions for "sloshing" of the magma chamber to occur.They found that the foam layer completely collapsed when oscillations were near the resonance frequency of the fluid layer and when sloshing lowviscosity fluids surrounding large bubbles.In nature, the collapse of a foam would potentially release a gas slug and produce a magmatic eruption, or result in magma overturn and a delayed eruption.

Magma chamber to pluton
Crystal-rich magmas are relevant to the final stage of a chamber as it reaches thermal equilibrium with the crust.Bons et al. (2015) developed a simple one-dimensional model using the finite-difference method to simulate the vertical profile of the evolving crystal fraction of an initially liquid magma chamber.Based on equations that are analogous to the inviscid Burgers' equation (Burgers, 1974), their model applied the so-called "traffic jam theory", which has also been used to explain the spontaneous formation of traffic jams on a motorway, where cars interact and slow down due to locally dense traffic.Their models suggest that self-organization of crystals will occur in a cooling magma reservoir due to gravitational sorting of floating or settling crystals.The distance the crystals travel depends on parameters such as melt viscosity and cooling rate.As the crystals interact barriers may form which instigate layers to develop, and each layer then undergoes similar crystal sorting.This could explain rhythmic layering but also larger-scale zonation observed in large igneous bodies.Bergantz et al. (2017) used a discrete-element computational fluid dynamics simulation to study the intrusion of basalt into an olivine-basalt mush by considering crystalscale force chains, crystal transport, and mixing processes.They found that the conditions for mechanical "lock up" were not uniquely identified by a specific crystal volume fraction.In comparison, Parmigiani et al. (2017) found that buoyancy-driven outgassing at the pore scale was most efficient at crystal volume fractions between 0.4 and 0.7, but that significant outgassing would be required to happen at higher crystal fractions and so require veining and/or capillary fracturing.
Examples in which large, high-viscosity magma bodies have been studied experimentally in the laboratory are relatively rare.Pluton emplacement has been modelled as a large body of viscous fluid which plastically deforms its surroundings, and experiments have focussed on pluton geometry and how space is accommodated in the lithosphere.Perhaps one of the earliest analogue experiments to study granitic pluton emplacement was carried out by Ramberg (1970), who used a combination of clay, putties, wax-oil mixtures, plates of concrete, and aqueous solutions to simulate diapiric ascent of fluid-like magmas through rock layers with differing competency (see Tables 1 and 2).In his experiments, Ramberg used a centrifuge model arrangement capable of reaching an acceleration of 4000 × g.Roman-Berdiel et al. (1997) and Roman-Berdiel (1999) studied granite emplacement by injecting low-viscosity Newtonian silicone putty into a tank of sand.Processes such as granite intrusion under the influence of tectonic stresses (Mazzarini et al., 2010), interaction with coincident faults and fractures during transpression (Benn et al., 1998), and strike-slip (Corti et al., 2005) motion have also been considered.Laboratory experiments by Roman and Jaupart (2017) document the post-emplacement behaviour of mafic intrusions as they spread at a density interface and founder as they become denser than their surroundings.They found that the dense magma body develops into different geometries depending on its aspect ratio, either creating a "teardrop" or "jellyfish" morphology with Rayleigh-Taylorlike instabilities at the intrusion margin.These experiments have collectively found that magma intrusion rate, mechanical heterogeneity of the crust, and tectonics play an important role in controlling the level of magma accumulation in the crust and its ability to reach the surface.

Testing magma chamber models
The numerical simulations and laboratory models of magma chamber dynamics can be tested against modern case studies with recent volcanic activity (e.g.Mt Pelée; Annen et al., 2008, see Fig. 7a) but also against the rock record to help interpret whether an exposed pluton was once a magma chamber or several discrete and small magma bodies (e.g. the Torres del Paine intrusive complex, Chile, Leuthold et al., 2012; and the Tuolumne Intrusive Suite, Sierra Nevada, USA, Coleman et al., 2004).Increasingly complex finite-element modelling techniques are being developed to account for the evolving thermo-mechanical and chemical processes associated with magma chamber recharge events.Such techniques are applied to volcanic centres that are experiencing periods of unrest to infer characteristics of the magma chamber including depth, overpressure, volume change, and shape (e.g.Hickey et al., 2016; see Fig. 7b).
One current limitation of numerical models of magma chambers is that they are predominantly static and often do not consider magma injection mechanisms; however there are exceptions.The dynamic numerical simulations of Longo et al. (2012) consider magma convection and mixing within a chamber that is injected by a buoyant magma.They relate the resulting pressure changes they model due to magma injection to the generation of ground displacements that could be detected at the surface.There is scope for more interaction between numerical and analogue modellers of magma chambers, with great potential to advance our understanding of this dynamic and yet highly enigmatic component of volcanic systems.

Sheet intrusions
Magma transport through the crust is facilitated by a series of interconnected magma-filled sheet intrusions called dykes and sills (see Fig. 1).Together these comprise a volcanic plumbing system that stores magma at depth but also can directly feed eruptions at the surface.A key assumption in many models is that sills are fed by dykes.The modelling approach depends on the assumptions on the controls of magma intrusion: (1) magma intrusion is modelled as a hydraulic fracture using the principles of linear elastic fracture mechanics and propagation is driven by fracturing of the host rock, or (2) magma intrudes as a viscous indenter and the growth dynamics are governed by the plastic deformation of the host and fluid properties of the magma.Scaling of magma intrusion models incorporates aspects of fluid flow but also host-rock deformation (Merle, 2015).The Reynolds number is typically modelled as low Re, meaning laminar flow.Host deformation is either scaled by the Hubbert number (the ratio of gravity stress to cohesion, and an important consideration in brittle materials) or the Ramberg number (the ratio of gravitational to viscous forces, and particularly relevant for scaling ductile materials).Geometric scaling compares, for example, the aspect ratio of the sheet intrusion (ratio of length to thickness) in the model and in nature.End-member host rheologies capture a wide spectrum of dynamic behaviour of sheet intrusion growth, from hydraulic fractures to viscous indenters.

Numerical models of sheet intrusions
A comprehensive review by Rivalta et al. (2015) explores the state of knowledge regarding dyke propagation models that use (1) Weertman theory (1971), in which the dyke is modelled as a buoyant magma-filled fracture; (2) lubrication theory (Spence et al., 1987;Lister, 1990), in which the dyke propagation is controlled by the flow of magma; or (3) a combination of both.A commonly used numerical approach to study dyke propagation is the boundary element method (BEM) (e.g.Dahm, 2000;Muller et al., 2001;Maccaferri et al., 2011), which considers the coupling between magma pressure and rock deformation, using analytical solutions for elementary dislocations to represent a pressurised, propagating crack.Dyke propagation in a stress field, controls on dyke trajectory, and tendency to form sills have been well studied numerically (e.g.Maccaferri et al., 2011;Barnett and Gudmundsson, 2014).Three-dimensional finite-element models are rare as they require re-meshing of the entire domain and so are more computationally demanding than the two-dimensional BEM approach, in which re-meshing requires only that new elements are added to the dyke tip (see Fig. 8a).The future of numerical model approaches to study dyke propagation is towards three-dimensional simulations to account for the complexities that are apparent in field geology studies and geophysical surveys.Numerical models of sills and laccoliths are often steady state (equilibrium, static, time invariant), two-dimensional, and axisymmetric.Malthe-Sørenssen et al. ( 2004) and Zhao et al. (2008) implement discrete-element models to study sill emplacement, accounting for fracturing of the host rock and considering sills and laccoliths of any size.The evolving intrusion geometry and surrounding deformation of the host material has more recently been studied by idealising deformation of the rock overburden as bending a stack of thin elastic plates, following the approach first outlined by Pollard and Johnson (1973), which is applicable to relatively shallow intrusions.Bunger and Cruden (2011) expanded the thin elastic plate theory to include fracture propagation criteria, fluid flow and the weight of the magma to explain the progression of the intrusion geometry from a bell-shaped geometry to flat top and steep-sided laccolith to thin disclike morphology of large mafic sills over time (Fig. 8b).In comparison, Michaut (2011) models shallow magma intrusions using non-dimensionalisation of the flow equation to describe magma spreading beneath an elastic crust, finding that the characteristic intrusion length depends on the elastic properties of the overburden and the characteristic intrusion thickness depends on the magma properties and the injection rate.The model of Galland and Scheibert (2013) accounts for axisymmetrical uplift both above and outside the intrusion, and it has recently been used to invert for laccolith dimensions and depth associated with surface deformation at Cordón Caulle volcano during a rhyolite eruption in 2011 (Castro et al., 2016).

Elastic deformation
Gelatine is a viscoelastic material that deforms almost ideally elastically at low temperatures and low concentrations (Kavanagh et al., 2013).It has been used as an analogue host material for magmatic sheet intrusions since Hubbert and Willis (1957)  (Fig. 3a).For decades since, this material has been used to investigate the dynamics of magma intrusions considering a large range of parameters such as density contrasts between magma and host rock (e.g.Fig. 9a; Taisne and Tait, 2011;Takada, 1990), mechanical layering of host rocks (e.g.Rivalta et al., 2005;Kavanagh et al., 2006), and impact of a stress field (e.g.Menand et al., 2010;Daniels and Menand, 2015).
The orientation of dykes in nature suggests they can be strongly influenced by the stress field in which they propagate.Fiske and Jackson (1972) used an unsupported gelatine solid to study the impact of gravitational forces on the trajectory of magma injections in the crust applied to volcanic rifts in Hawaii (Fig. 3b).When water was injected, the dyke that formed moved laterally following the ridge axis and was oriented perpendicular to the least compressive stress direction.Subsequent gelatine studies have investigated dyke injection beneath or into a conical edifice (e.g.Hyndman and Alt, 1987;McGuire and Pullen, 1989;Muller et al., 2001;Kervyn et al., 2009; see Fig. 9c), the impact of the formation of a collapse scarp on dyke orientation (e.g.Walter and Troll, 2003), and dyke injection into gelatine under extension (Daniels and Menand, 2015) and compression (Menand et al., 2010).In all cases, the experimental sheet intrusions' propagation path was modified when entering a new stress field, realigning with the new σ 1 direction.
Due to the photoelastic properties of gelatine, stress in the host material associated with the formation and growth of fluid-filled dykes and sills is shown in polarised light (e.g.Kavanagh et al., 2017; see Figs. 9a and 10a), and the evolution of incremental and finite strain in the deforming gelatine due to dyke, sill, and hybrid intrusions has been quantified using DIC (Kavanagh et al., 2015(Kavanagh et al., , 2017; see Fig. 10b, c).Mechanical heterogeneities such as layering (e.g.Rivalta et al., 2005;Kavanagh et al., 2006) and the presence of discontinuities such as faults or joints (Le Corvec et al., 2013) has been shown to influence the propagation of a fluid-filled crack in an elastic gelatine host material and cause the transition of a dyke into a sill.The results show there are significant (up to 60 %) decreases in strain and stress around a feeder dyke as a sill forms in an elastic material (Kavanagh et al., 2015(Kavanagh et al., , 2017; see Fig. 10c), and this can be directly linked to a decrease in fluid pressure when the sill forms (Kavanagh et al., 2015(Kavanagh et al., , 2017)).Layered gelatine experiments have also been used to study laccolith emplacement whereby viscous grease was injected into layered gelatine with a lubricated interface (Pollard, 1973;Pollard and Johnson, 1973;Johnson and Pollard, 1973).

Compacted granular materials
Compacted finely grained silica flour has been used as a hostrock analogue in which magma intrusion is modelled as a viscous indenter as the host deforms plastically (e.g.Abdelmalak et al., 2012;Mathieu et al., 2008;Guldstrand et al., 2017;Schmiedel et al., 2017).This material can fail both in tension and in shear due to its non-negligible cohesion (Abdelmalak et al., 2016).Other comparable granular materials that have been used to study sheet intrusions include ignimbrite (Mathieu et al., 2008), diatomite (Gressier et al., 2010), and dry plaster powder (Kervyn et al., 2009), which, in general, are finely grained frictional, cohesive, and variably permeable (see Table 1 for material properties and Table 2 for host rock and magma analogue combinations).Factors such as the surface deformation associated with sheet intrusion (Galland et al., 2016;Guldstrand et al., 2017), mechanical strength of the host rock (Schmiedel et al., 2017), mechanical layering (Gressier et al., 2010), and ambient stress field (Galland et al., 2006) have been considered using these materials.As they are opaque this requires solidifying fluids to be injected (often Vegetaline; see Table 2) and then either the resulting intrusion is excavated postemplacement and its dimensions linked to surface deformation (see Fig. 10d for an example experiment set-up; Galland et al., 2014), or a thin quasi-two-dimensional tank is used to show a cross section through the experiment as intrusions form (e.g.Abdelmalak et al., 2012).A variety of laboratory intrusions have been created in compacted silica flour experiments, ranging from cone sheets to dykes and sills (Galland et al., 2014; see Fig. 9d).Galland et al. (2016) used structure from motion to reconstruct the three-dimensional geometry of the excavated solidified intrusions and geo-referenced these to the surface deformation that was produced.For sill emplacement, Galland (2012) noted the symmetrical updoming of the initially flat experiment surface in response to the sill emplacement, with the complexity of the subsurface intrusion that was later excavated being comparable to the complexity of the surface deformation it caused.The formation and growth of visible fractures, lateral and vertical displacement, and calculation of shear strain based on surface changes with time was mapped using DIC.

Solidification and viscosity effects
Experimental work has shown that the viscosity of the fluid within sheet intrusions can impact the geometry of dykes and sills as they propagate in their host material.Table 1 lists various magma analogues that have been used in dyke laboratory experiments, including viscous liquids such as golden syrup and honey and solidifying liquids such as vegetable oil (Vegetaline).Solidification within dykes has been studied in gelatine experiments (e.g.Taisne and Tait, 2011; Fig. 9b) and shows a transition in propagation behaviour of the dyke tip from continuous to step-wise, with progressive stalling, infla-tion, and then breakout due to breach of an insulated and relatively low-viscosity fluid from the intrusion interior through a cooled margin.Chanceaux andMenand (2014, 2016) formed sills between gelatine layers by injecting solidifying Vegetaline directly into an interface (Fig. 10e-f).They found that the sill tip region has the tendency to become segmented, especially when the injected fluid is viscous or has become more viscous due to solidification.

Testing magma intrusion models
There are several challenges that mean testing magma transport models is not straightforward.By their very nature, magma intrusions are subsurface features and so cannot be directly observed during their formation or when they are active.Insight into active intrusion processes in nature is typically interpreted based on analysis of surface deformation thought to be related to magma movement and in combination with seismic data that are inferred to result from intrusion-related rock fracturing (e.g.Sigmundsson et al., 2014).Due to these limitations, there is much discussion in the literature regarding how to model dykes and how horizontal sheet intrusions (sills) relate to the construction of larger igneous bodies such as laccoliths and magma chambers.It is likely that magma viscosity (and composition) will affect the mechanisms of sheet intrusion, with highly viscous magmas more likely to cause shear failure of their host material than tensile opening.It is also likely that a spectrum of host rheology exists in space and time during sheet intrusion emplacement, and so future work should involve exploring complex materials that span the end-member host rheologies.Numerical models of a penny-shaped hydrofracture propagating in an infinite elastic material (e.g.Savitski and Detourney, 2002) have been used to interpret the results of hydraulic fractures in gelatine experiments.Using measurements from gelatine experiments (e.g.intrusion dimensions, injection flux, fluid viscosity, host material Young's modulus, and Poisson's ratio), it has been shown that dyke and sill propagation occur in a toughness-dominated regime, in which the fracture properties of the host material control the dynamics (e.g.Menand and Tait, 2002;Kavanagh et al., 2017).However, in some cases sill growth in gelatine experiments has been better explained by viscosity-controlled dynamics (e.g.Kavanagh et al., 2006;Chanceaux and Menand, 2016).
There is great potential to use laboratory models and numerical models of magma intrusion in combination to assist the development of inversion methods to characterise magma intrusion geometry and depth in nature.The numerical models that calculate intrusion geometry and depth using ground deformation measurements such as GPS and InSAR at active volcanoes (e.g.Fukushima et al., 2005) could be tested on laboratory models of magma intrusion in which the volume, depth, and geometry of the experimental dyke or sill is known.Considering additional complexities such as natural topography, tackling non-axisymmetric intrusion geometries, and assessing the impact of inelastic deformation of the host (e.g.Scheibert et al., 2017; see Fig. 8c), heat exchange with the surroundings (e.g.Thorey and Michaut, 2016), magma cooling and solidification, mechanical layering of the host-rock, local or regional stress perturbations, the presence of weak discontinuities such as faults or fractures, and pressure variations within the intrusion all result in a non-unique set of best-fit simulations when applying models, and so constraining which parameters to include and exclude needs to be evaluated on a case-by-case basis.Comparison with data from experiments that monitor surface deformation (e.g.Kavanagh et al., 2015;Galland, 2012;Guldstrand et al., 2017;Tortini et al., 2014) will enable the validation and improvement of inversion models, with the identification of those experimental parameters that are modelled well and those which are not to help guide future research for hazard assessment at active volcanoes.

Lava lakes
Lava lakes are effusions of lava either at the top of an open conduit (e.g.Halema'uma'u lava lake, Kīlauea), or a ponded area of an active lava flow (e.g.Kīlauea Iki 1959-1960, Kīlauea, Hawaii;Richter et al., 1970).There are several active lava lakes in the world, each with a range of hazards, from gas plumes and outpourings of lava to explosions which could occur without warning.The dimensionless numbers used to study lava lakes and their connected conduit are the Prandtl, Rayleigh, and Reynolds numbers (e.g.Harris, 2008).Both low and high Re are explored, but the transport number Te which is dependent on Re is also explored (e.g.Huppert and Hallworth, 2007).

Lava lake surface morphology
Techniques such as time-lapse photography have shown that lava lake surfaces in nature are highly dynamic parts of a volcanic system (e.g.Orr and Rea, 2012), and spreading of the gradually cooling lava lake surface has been linked to convection in the underlying lake and conduit.Karlstrom and Manga (2006) used molten paraffin wax to study lava lake dynamics, monitoring the surface of the wax with infrared cameras.The wax surface was cooled before a partially submerged bar pulled the crust apart at a constant velocity along an incision in the wax surface.This formed zigzag rifting morphologies reminiscent of structures described at natural lava lakes and enabled the calculation of the spreading rate, crustal thickness, and yield stress of the crust and thus the strength of the convective forces acting upon the underlying magma.Harris (2008) used numerical modelling to study convection of a molten lava lake fed by a conduit.In his models, an upwelling injection of hot, degassing, buoyant, and less-viscous magma rises through the conduit to the lava lake.Radiative heat loss and surface spreading then induced cooling and an increase in fluid density that caused downwelling of the magma, with the highest viscosity magma flowing back down the conduit.This model is supported by water and/or glycerine models of Huppert and Hallworth (2007), who showed that different density fluids may flow past one another during exchange flow in the conduit.Beckett et al. (2011) developed this further using syrup models exploring the impact of viscosity ratios between the different density fluids on the exchange flow.Molina et al. (2012) studied the lava lake at Mt Erebus using a numerical model of a 4-10 m diameter conduit filled with melt plus crystals.They found that larger conduits allow for greater convection rates and maintain a high enough temperature to remain viscous for extended periods of time (30 years).From their model, the surface spreading rate was calculated at between 10 −6 and 10 −5 m s −1 ; however this modelled value is dramatically slower than what is observed at Erebus, suggesting that a simple melt-plus-crystal magma does not capture the full dynamics and that a gas phase should be accounted for in future models (Molina et al., 2012).

Lava lake surface level variations
Fluctuations in lava lake level have been associated with gas exsolution.Witham et al. (2006) carried out a series of experiments that released air from a deep "chamber" into a cuboidal conduit (1 × 1 × 18 cm) of water attached to an approximately cubic (14.1×14.1×15cm) "lava lake".Gas was released into the base of the conduit using a compressor; this decreased the water-air density, causing the bubbly mixture to rise into the surface reservoir due to buoyancy, resulting in an increased lava lake level.Gas was then released from the water at the surface of the higher reservoir, progressively increasing the hydrostatic pressure.Eventually the hydrostatic pressure of degassed water in the lake exceeded the pressure from below, preventing further rise of gas-rich water and resulting in collapse of the conduit, fluid flow back down into the chamber, and lowering of the lake level.These analogue experiment results show that rising lava lake levels in nature could be explained by periods of increased gas emission from the chamber through the conduit, and that decreases in lake level could occur when the magma static pressure in the overlying magma column exceeds the pressure of the rising magma.

Lava domes
Lava domes are effusions of degassed, highly viscous, silicarich magma that accumulate at volcanic vents.Their emplacement can cause the build-up of gas and pressure in the conduit, increasing the potential for explosive eruptions or the formation of pyroclastic density currents.Modelling lava www.solid-earth.net/9/531/2018/Solid Earth, 9, 531-571, 2018 dome emplacement and stability is key for identifying thresholds for collapse and therefore for assessing the potential risk of such events.Aspects of lava dome emplacement that have been studied in laboratory experiments include morphological variations due to topography (e.g.Lyman et al., 2004), magma rheology (e.g.Griffiths and Fink, 1993;Fink and Bridges, 1995), and the preservation of flow fabrics using magnetic fabrics (Závada et al., 2009).Lava dome models have been scaled by considering the Reynolds number (Re) ensuring laminar flow and the Bingham number (B), which quantifies stresses within the dome (Balmforth et al., 2000).

Dome morphology
Laboratory models of lava domes have largely focussed on the influence of lava rheology on dome morphology.Griffiths and Fink (1993) investigated the progressive spreading of lava domes by effusing liquid PEG 600 wax into a tank of cold sugar solution with a horizontal base.The temperature gradient between the wax and solution caused the onset of solidification, and the lava viscosity had a large influence on the morphology of the dome that was formed.Fink and Bridges (1995) found that pulsating the wax effusion and decreasing its temperature resulted in predominantly vertical growth of the dome rather than flow away from the vent, and so the length of lava domes could be explained primarily by variations in effusion rate.Several lines of evidence suggest that lava dome rheology in nature is more complex than a simple temperaturedependent Newtonian fluid.Balmforth et al. (2000) carried out numerical simulations of lava dome growth and evolution using a Herschel-Bulkley rheology, which were compared with kaolin-water slurry models.They found that the yield stress acting in the dome is important in determining dome morphology; however the combined effects of shear thinning and yield stresses were difficult to distinguish.Experimentally, Blake (1990) used kaolin-water slurries to model dome morphology in a Bingham fluid to identify yield stresses occurring within domes effused onto horizontal planes.Their domes were parabolic in nature with surface lineations that spiral out from the centre and are inferred to be areas of shear localisation, which allows for stable dome growth.Griffiths and Fink (1997) used a PEG-kaolin mixture to study lava dome morphology with the kaolin powder converting the fluid from a temperature-dependent Newtonian fluid to a Bingham fluid.These analogue lava domes produced spines and irregular breakouts of wax (Fig. 11a) due to the yield strength of the magma analogue.Lyman et al. (2004) used a similar mixture to investigate the impact of slope and effusion rate on dome morphology.They found that contrasting dome morphologies (e.g.platy, spines, lobes) were associated with extrusion onto a surface at different slope angles, but that effusion rate had the greatest impact on dome morphology.These results can be compared to lava dome morphologies in nature, such as Wilson Butte in California (Ly-man et al., 2004), to calculate the effusion rate of prehistoric domes.
Numerical models of dome growth have explored the impact of lava rheology and spatially varying mechanical properties on lava dome morphology, the development of spines, and the onset of dome failure.Some of the early theoretical models of lava domes were quasi-static and developed by Iverson (1990), who conceptualised lava domes as pressurised magma held within a brittle shell.He found that the thickness and tensile strength of the shell and the excess pressure of the enclosed magma and gas could be combined to give a dimensionless number D that governs the dome shape.Hale and Wadge (2003) built upon this to construct a finiteelement model that includes the rheological evolution of the magma due to degassing.They found that the morphology of modelled domes fits better with natural examples when an evolving rheology is considered.Hale (2008) then included an independently evolving talus slope of the dome, later validating the model by simulating the development of the dome at Soufrière Hills Volcano, Montserrat (Hale et al., 2009).More recently Husain et al. (2014) have used finite-element modelling to investigate the effects of extrusion rate and rheology on the growth of lava domes, determining that an evolution in the rheology allowed the softer more ductile magma in the dome core to become stiffer, enabling extrusion from the dome as spines.They also found a strong relationship between extrusion rate and the mechanical properties of the lava which impacted lava rheology.

Internal deformation of lava domes
Internal flow patterns within lava domes in nature have been inferred from crystalline and bubble fabrics within the crystalline lava, thus providing insight into the processes occurring within a dome during formation.Závada et al. (2009) studied magnetic fabric development within lava domes by effusing plaster of Paris seeded with magnetite particles from a point source, injecting with increasing pressure onto a deformable surface of sand (Fig. 11b).The plaster of Paris and magnetite mixture behaves as a shear thinning fluid and was allowed to solidify once extruded.The solidified dome was then cut into slices and oriented samples drilled for analysis by applying anisotropy of magnetic susceptibility (AMS) to quantify the direction and intensity of any fabric that was developed by the magnetite particles during the extrusion of the lava dome.They found more concentrated suspensions with higher viscosity created complex dome structures that had relatively steep sides, akin to lava domes commonly observed in nature.
Post emplacement cooling and alteration of the dome from hydrothermal fluids was modelled by Ball et al. (2015), who simulated the effect of heat and water flow though cooled lava domes over prolonged timescales (100 years).They determined that alteration is most likely to occur at boundaries between different parts of the dome, and that temperature greatly affects the potential for alteration with faster cooling rates reducing the likelihood of deep alteration within the dome.

Lava flows
Laboratory and numerical modelling has been extensively applied to investigate lava flow emplacement (Fig. 1).Scaling of lava flow models considers several aspects, such as the mechanism of heat loss via the Peclet number (Pe) and dimensionless timescales ( ) for solidification (e.g.Fink andGriffiths, 1990, 1998).
In the laboratory, the emplacement dynamics and morphology of lava flows has been investigated using Newtonian fluids (e.g.glucose syrup; Stasiuk et al., 1993) and more complex fluids that account for cooling and crystallisation and develop a solidified crust during flow (e.g.PEG wax; Hallworth et al., 1987;Fink and Griffiths, 1990;Gregg and Fink, 1995).These experiments model variations in heat flux, thermal gradients, and cooling on the temporal and spatial variation in lava flow viscosity, extrapolating on the impact these factors have on run-out length and flow morphology, for example.Lavas have also been modelled in the laboratory as a particle suspension, with experiments showing that increasing particle volume fraction (Soule and Cashman, 2005;Castruccio et al., 2014) and particle size (Del Gaudio et al., 2013) increases lava viscosity and can affect lava flow morphology.High-concentration particle suspensions produce low flow velocities, shear localisation, and subsequent break-up of the flow surface, causing transition from pahoehoe-like to aa-like morphologies that are reminiscent of natural flows (Soule and Cashman, 2005).

Lava flow dynamics
The cooling mechanism of lava flows impacts the morphology and run-out length of the flow and therefore the hazards it may pose to surrounding populations.Recent work by Garel et al. (2012Garel et al. ( , 2014) ) has studied, in detail, the effects of cooling on the morphology of flows in the laboratory and numerically.They effused silicone oil (Garel et al., 2012) and PEG P3515 (Garel et al., 2014) onto a horizontal piece of polystyrene, which insulates the base of the flow, and developed a numerical model that considers both surface and basal cooling by scaling thermal boundary conditions and radiated power.They showed that cooling is primarily controlled by the effusion rate of the fluid, and that it is possible to estimate the effusion rate from observation of the surface thermal signal.Their results correlate with natural examples of effusion rates; however there is a large amount of uncertainty attributed to the complexities associated with natural lava flows (Garel et al., 2012(Garel et al., , 2014)).

Lava flow levees and solidification
Critical in the cooling of lava fields is the development of lava levees, crust formation, and progressive breakout.This has been investigated in the laboratory using paraffin wax, where the progressive cooling of the hot, liquid wax causes levees to form and channelisation of the flow (Blake and Bruno, 2000;Miyamoto et al., 2001;Nolan, 2014).Crust formation over the cooling flow surface insulates the molten wax and creates tube-fed flows, and blockages or restrictions in the tube-fed flow of wax to the flow front lead to flow inflation and eventually breakout from the crust.Blake and Bruno (2000) used PEG wax experiments to demonstrate the link between lava effusion rate, lava viscosity, and strength of the chilled crust, which impacts how and where breakouts from lobate structures occur.Karlstrom and Manga (2006) used spreading paraffin wax experiments to study the morphological transition from pahoehoe to aa flows due to breakouts from the cooled, spreading crust (also see Sect.7.1 on lava lakes).
Solidification and development of columnar jointing in lava flows has been modelled using corn starch slurries that are placed under heat lamps to allow the water to evaporate away (Goehring and Morris, 2005;Müller, 1998).The loss of water was used as an analogue for heat loss within lavas; as the starch dries out it shrinks, resulting in cracks forming and propagating through the material (Fig. 12b).The morphology of the vertical columns formed within the analogue lava correlates well with the morphology of columns in natural lava lakes and ponded lava flows such as in Hawaii (Goehring et al., 2006;Müller, 1998) or the Giant's Causeway in Northern Ireland (Goehring and Morris, 2005).However, further rhewww.solid-earth.net/9/531/2018/Solid Earth, 9, 531-571, 2018 ological studies to better understand the material properties are needed to improve scaling these experiments to nature.

Substrate erosion
Field observations of erosion channels within lava tubes suggest that assimilation of the lava substrate can occur when lava flows are emplaced with high heat flux or flow over substrate with a low melting temperature.Huppert and Sparks (1985) investigated the development of thermal erosion channels in komatiite lava flows by pouring hot water onto a slab of PEG 1000 wax; Komatiite lavas are thought to have had unusually high heat flux, and so thermal erosion of their substrate is expected to have been an important process in the development of these ancient flows.Kerr ( 2001) used theoretical models alongside molten PEG 600 effused onto an inclined sheet of solid PEG 600 (Fig. 12a) to investigate how the thermal profile of lava flows evolves both spatially and temporally.His experimental results agreed with the theoretical models, which showed that there is a critical thickness range at which chilled margin formation at the base of the flow ceases and erosion begins, depending on the initial temperature of the lava; for basaltic lavas on Hawai'i, this range is 7.3 to 34 cm after a period of 0.21 to 4.6 days.

Flow indicators in lavas
When studying ancient solidified flows in nature, crystal distribution and stretched bubbles have been used to infer flow direction.The preservation of flow indicators in solidified lavas was investigated in analogue experiments using layered viscous silicone to model internal strain within extruding and spreading fluids (Gilbert and Merle, 1987).The experiments showed that in channelised flows, or at the base of a lobe, the lava flow trajectory indicators could be both parallel and perpendicular to each other in the upper portion of lobes.When applied to lava flows in nature, such observations can explain emplacement mechanisms and possibly account for variations in deformed bubble and crystal shapepreferred orientations compared to AMS fabrics in different parts of the flow (e.g.Caballero-Miranda et al., 2016).

Application of numerical lava flow models
Numerical models of lava flow inundation and behaviour have a range of complexities, from deterministic or computational fluid dynamics models designed to replicate emplacement behaviour to the more statistical stochastic models designed to estimate the probability of inundation (see Costa and Macedonio, 2005;Cordonnier et al., 2016;Dietterich et al., 2017).Lava flow numerical models have been used to simulate edifice growth from the accumulation of multiple lava flows (Annen et al., 2001), the insulating properties of lava tubes (Keszthelyi, 1995), cooling of pahoehoe lavas (Keszthelyi and Denlinger, 1996), and formation of lava levees (Quareni et al., 2004), and from a hazard perspective, flow length, rates of emplacement, and inundation areas.For hazard and risk assessment purposes, simulated lava flows are emplaced over a digital elevation model to predict the inundation pathways of flows.Applied models range from complex 3-D models, to simplified 2-D models (e.g.Costa and Macedonia, 2005), cellular automata models (e.g.Connor et al., 2012;Rongo et al., 2016), simplified 1-D models (e.g.Harris et al., 2016), and stochastic models (e.g.DOWN-FLOW; Tarquini and Favalli, 2011), with different model types appropriate for different eruptive scenarios (Dietterich et al., 2017).As is the case for other geophysical flows, 2-D models are based on shallow-water equations describing the conservation of depth-averaged velocity, temperature, and flow thickness.Examples of such models include VolcFlow (Kelfoun and Vargas, 2015) and RHEOLEF (Bernabeu et al., 2016).Cellular automata models are commonly used in probabilistic hazard assessment of lava flow inundation (e.g.Miyamoto and Sasaki, 1997;Crisci et al., 2004;Favalli et al., 2005;Vicari et al., 2007;Connor et al., 2012;Rongo et al., 2016) (Fig. 12c).Such methods account for mass conservation, radiative heat exchange with the atmosphere, viscosity variation with temperature, and non-Newtonian flow behaviour.In these models, the area of interest is discretised into a regular grid of cells or points (in the example of SCIARA, cells are hexagonal in shape; Rongo et al., 2016), which can be inundated depending on a cell's elevation in relation to neighbouring cells.In Connor et al. (2012), the probability of a cell being inundated is dependent on the relation between the elevation of the empty cell and the thickness of the lava in neighbouring cells.Such models assume a given volume is erupted, with this volume being distributed amongst the inundated cells.An example of a simplified 1-D model is FLOWGO (Harris et al., 2016), which tracks the thermo-rheological evolution of a lava flow as it flows within a channel following the steepest line of descent.The model assumes Newtonian flow, modified for a Bingham fluid, and using a series of heat loss equations is able to calculate the cooling and crystallisation and as a consequence, flow viscosity and emplacement velocity.However, application of the model requires assumptions regarding the channel dimensions at vent to match effusion rate inputs, limiting interpretation of results.While the simplicity of the modelling techniques above means that they can be quickly and easily employed for rapid hazard assessment in the build-up or during an eruption, it is argued that this simplicity means that it is not possible for them to capture the distribution of lava with time (Bernabeu et al., 2016;Dietterich et al., 2017).Such behaviour requires the application of more complex 3-D modelling techniques, and increasingly, computational fluid dynamics tools such as OpenFOAM and FLOW-3D are being employed to simulate the complex processes associated with emplacement of volcanic flows (see Cordonnier et al., 2016 andDietterich et al., 2017).Both OpenFOAM and FLOW-3D are able to compute heat transfer and can simulate two-phase and viscous flows and a range of rheologies.These fully thermodynamical rheological models have higher computation requirements than the simpler examples but are considered more appropriate for the simulation of long-lived and cooling-limited flows (Dietterich et al., 2017).
10 Particle-laden flows Particulate-laden flows are one of the most hazardous phenomena associated with explosive volcanic eruptions with the potential to impact areas hundreds of kilometres from the source.Our understanding of pyroclastic density current and lahar processes is particularly reliant on numerical and laboratory modelling due to the hazardous nature of the phenomenon.

Pyroclastic density currents
The term pyroclastic density current encompasses a wide range of flows from dilute surges to dense flows, block-andash flows, and pumice flows (Branney and Kokelaar, 2002;Sulpizio and Dellino, 2008).Such flows occur as a mixture of particles and gas that is emplaced as a gravity current which propagates down the slopes of a volcano and is related to, for example, collapse of a volcanic column or lava dome.For the purposes of this overview, it is sufficient to consider a pyroclastic flow as divided into two main parts: a dense basal portion where movement is controlled by particle-particle interaction, and an overlying turbulent dilute region that is composed of ash and gas.Friction dominates at the basal portion of the flow, where the Savage number Sa is used to describe the collisional stresses related to frictional stress.The Darcy number Da describes flow fluidisation, whereby the effective permeability of the particle mixture decreases, enabling high pore pressures to develop, being greater than 1 if the region is fluidised.Above the basal layer of the flow, the particle concentrations are slightly lower and momentum is distributed by particle momentum (described by the collisional stokes number St c and the Bagnold Ba number, which is the ratio of collisional stresses relative to viscous stress).Within the dilute, turbulent portion of the flow, particle motion is influenced primarily by particle-gas drag as described by the Stokes number St. Other dimensionless parameters used in describing flow propagation are the Peclet and Reynolds numbers, the Froude number Fr (the flow speed relative to wave speed), and the dense-dilute transition number D D (the timescale of collisions relative to particle drag response time; Dufek, 2016).
It is not well understood how the portions of the pyroclastic density current interact; this complexity, and additional numerical challenges, have resulted in laboratory and numerical modelling techniques typically representing either the lower dense portion of the flow (e.g. the flow simulation software TITAN2D) or the dilute portion of the flow (e.g.Bursik and Woods, 1996;Andrews and Manga, 2011).

Laboratory investigations of pyroclastic density current dynamics
Some of the first laboratory experiments to investigate controls on pyroclastic density current propagation involved the injection of dense, sometimes particle laden, fluid into less dense fluid (Carey et al., 1988;Huppert et al., 1986;Sparks et al., 1993;Woods and Bursik 1994;Woods and Caulfield, 1992).In the example of Carey et al. (1988), buoyant plumes were produced by injection of particle-laden freshwater into saline water.Flows formed when the particle concentration of the injected fluid was large such that the density difference with the ambient fluid was reversed, leading to collapse of the plume.Woods and Bursik (1994) specifically focussed on the movement of flows on slopes and their interaction with topographic barriers (Fig. 13a).In this example, a dense fluid composed of a mixture of methanol and ethylene glycol was injected or released into water.Such experiments provided information on the control of entrainment on flow density, and controls of topography on both entrainment and sedimentation from a flow.More recently laboratory experiments have considered the effect of interstitial pore pressure on flow motion (Fig. 13b; Roche, 2012;Rowley et al., 2014) and the controls on co-ignimbrite plume formation from dilute flows using talc to represent finely grained ash particles (Andrews andManga, 2011, 2012).A key issue with regards to modelling volcanic phenomena is scaling, and over the past 5-10 years, significant effort has focussed on the development of so-called "large-scale" experimental set-ups to overcome this (Fig. 13c) (Dellino et al., 2007;Lube et al., 2015;Valentine et al., 2015) and try to more accurately reproduce observed dynamics.Such modelling also allows the use of natural eruptive pyroclastic materials to more accurately reproduce the physical dynamics, states, and relations that occur within real flows.

Numerical simulations of pyroclastic density currents
Numerical models that describe either the dilute or dense endmember flow are typically depth averaged and solve equations for conservation of mass, momentum, and thermal energy (e.g.Bursik and Woods, 1996).Such models are also steady state, and therefore do not account for changes in flow behaviour with time.An example of a numerical model that has been developed to account for both the tur-bulent upper layer and the dense layer is the transient model of Doyle et al. (2008).In this model, the dilute current is described by depth averaged, isothermal, continuum conservation equations, while the basal flow is modelled as a granular avalanche of constant density.
In addition to the production of separate models to account for different physical processes, numerical models of varying complexity also exist.Numerical models have been developed that consider dilute particle-laden gravity currents (e.g.Bursik and Woods, 1996;Dade and Huppert, 1996) to calculate properties such as velocity, temperature, and density of the flows.These calculations have a small number of parameters and as such involve a number of simplifications, allowing parametric studies to be conducted to understand the control of inputs on the modelled outputs.
The advancement of computational efficiency has enabled the development of supercomputer calculations, solving full Navier-Stokes equations for flows and fountains (Valentine et al., 1992;Neri and Dobran, 1994;Esposti Ongaro et al., 2012).These are elaborate computer codes, incorporating a large number of input parameters, and solving for a number of different phases, for example fluid (magmatic and atmospheric gases) and particles of different sizes and density simultaneously (Fig. 13d).While numerical studies of pyroclastic density currents still largely follow these two strands, open-source computational fluid dynamics programs (e.g MFiX and OpenFOAM) are being increasingly used to capture the multiphase behaviour that occurs within these volcanic flows.A comprehensive overview of the numerical models describing pyroclastic density currents is provided in Dufek (2016).

Lahars
Primary lahars form when an eruption causes melting of ice overlying the volcano, and eruptive products mix with meltwater to produce high-density mixtures of water and debris.Secondary lahars form as debris emitted by an eruption is mobilised after deposition, usually in relation to heavy rainfall between eruptive events.Similar to pyroclastic density currents, lahars can be described by two endmembers: stream flows and debris flows, but with an intermediate flow type called hyper-concentrated flows.Flow type can vary within an individual event both spatially and temporally, in association with changes in channel and underlying topography (e.g.Manville et al., 2013).Definition of lahar flow type is dependent on the concentration of particles, with stream flows representing those with low particle concentrations and debris flows with high particle concentration.The flow endmembers have very different rheology, and as a result the application of both laboratory and numerical modelling techniques is affected by similar challenges to those for pyroclastic density currents.Given similarities in flow behaviour, a number of the dimensionless parameters used to describe pyroclas-tic flows are also common to lahars, in particular relating to basal fluid pressure.

Laboratory modelling of lahars
Lahar laboratory experiments have largely considered debris flows, i.e. those with high particle concentration.As for pyroclastic density currents, both small-and large-scale (metre and tens of metres, respectively) laboratory experiments have been conducted.Iverson (2015) showed that experiments of debris flows are particularly susceptible to scaling issues and as a result, the United States Geological Survey has developed a 95 m long flume that allows the release and flow of up to 10 m 3 of water-saturated sediment onto a bed with variable roughness to closely mimic conditions in natural debris flows (Iverson et al., 2010(Iverson et al., , 2011)).Such experiments are used to show the relation between volume of flow and run-out and investigate the process of bulking whereby substrate is incorporated into the moving flow and changes in flow velocity and behaviour with variation in slope.

Numerical simulation of lahars
There are relatively few numerical models that simulate lahars.Numerical modelling of lahar dynamics is non-trivial due to spatial and temporal variation in flow behaviour, for example rheology (see Manville et al., 2013).Perhaps the most commonly used numerical model to simulate lahar emplacement is LAHARZ (Schilling, 2014), a computational model that uses empirical relations of past inundation events to forecast inundation for a given future event.More complex numerical models applied to lahars include a version of TITAN2D that accounts for both particles and fluids (Pitman et al., 2003;Williams et al., 2008), the commercial hydraulic model Delft3D (Carrivick et al., 2009), GIS-based models (Darnell et al., 2012), and application of models more typically applied to water floods, e.g.LISFLOOD.No single model can account for the all the different phenomena described above, due to the complex interactions between the solid and liquid phases.As such, research into the rheology of lahars is required to provide a description of the underlying physics to be utilised in numerical models.

Volcanic plumes
Much of our current understanding of plume dynamics can be traced back to the 1950s and the work of Morton et al. (1956) (see Fig. 3c).While not focussing specifically on volcanic plumes, their study has been the basis for much of the subsequent research into volcanic plumes using laboratory modelling and numerical models to reproduce behaviour observed in nature.Key dimensionless parameters for volcanic plumes are the Reynolds number Re and the Richardson number Ri, buoyancy and stratification parameters F o and G, respectively, and the Brunt-Väisälä frequency N.

Laboratory modelling of volcanic plumes
Traditional laboratory modelling of volcanic plumes has involved the injection of a less-dense fluid (e.g.freshwater or methanol) into a tank of denser fluid (e.g.saline fluid or ethylene glycol, respectively).The difference in densities between the two fluids allows the less-dense fluid to rise through the dense fluid, reproducing those characteristics associated with buoyant plume rise (see Fig. 14a; Carey et al., 1988).Injection into a stratified fluid enables not only modelling of plume rise but also the dynamics of plume spreading once the injected mixture reaches neutral density (e.g.Carey et al., 1988).Variation in the injection rate provides firstorder information on plume dynamics, and in particular on the controls on column collapse (Woods and Caulfield, 1992;Kaminski et al., 2005).In some examples, the injected fluid is particle laden to investigate the additional effect of particle sedimentation and re-entrainment on buoyant columns (Ernst et al., 1994;Carey et al., 1988), such that experimental observations can be related to those seen in ash fallout deposits in the field.In particular, Ernst et al. (1994) looked at the effect of wind on volcanic plumes, reproducing the bifurcation of the plume.These experiments are based around buoyant theory, and therefore overlook the processes occurring in the jet region of the plume, where dynamics are controlled by the upward velocity of material as it is ejected from the vent.Such laboratory experiments underlie the onedimensional numerical models currently used to investigate volcanic plume behaviour and inform inputs for ash-dispersal models widely used today.
A key issue with laboratory modelling of volcanic plume dynamics is scaling, particularly when considering that plumes in nature are injected into a stratified and turbulent atmosphere and are affected by local weather patterns.To account for such scaling issues, large-scale experiments are increasingly applied to investigate eruption processes (e.g.Dellino et al., 2014;Fig. 14b).These experiments, which often use natural materials, have enabled analysis of the effect of vent conditions and processes on subsequent eruption behaviour, and in particular have been used to estimate the rate at which air is entrained into the rising plume (the entrainment coefficient), which is a key input parameter for numerical modelling of plume rise (Costa et al., 2016).

Application of numerical models for simulating plume rise
Numerical models of volcanic plumes serve two main purposes: (1) to provide input information (for example plume height and mass flux of ash into the atmosphere) for ashdispersal models, and (2) to investigate the controls on these parameters.Two types of numerical model are available: onedimensional (integral) models and multicomponent multiphase three-dimensional models.

One-dimensional plume models
One-dimensional or integral models (Fig. 14c) are commonly used for defining source parameters for ash-dispersal models, largely because they are computationally inexpensive and results can be acquired quickly.These models (e.g.Bursik, 2001;Mastin, 2007;Barsotti et al., 2008;Folch et al., 2016;de' Michieli Vitturi et al., 2015) account for conservation of mass, momentum, and energy and are largely based around the model developed by Morton et al. (1956), with the constituent equations modified for application to the volcanic example by Wilson and Walker (1987), Sparks (1986), andWoods (1988).Since the 2010 eruption of Eyjafjallajökull, there has been increased emphasis on development of models that are able to account for the effects of wind on a rising plume (e.g.Woodhouse et al., 2013;Degruyter and Bonadonna, 2013).The models assume that the emitted gas and particles are in dynamic and thermal equilibrium; an approximation that is reasonable for dilute, finely grained plumes but is less appropriate in the jet part of the plume.Entrainment coefficients are key inputs for plume models, and these parameters have been the focus of much research in recent years.Entrainment in one-dimensional models is captured using two additive entrainment parame-ters, one accounting for radial entrainment associated with the incorporation of ambient air by turbulent eddies at the plume edge, and the second accounting for the effect of wind on air entrainment.The first coefficient has been well defined using observations from laboratory experiments of turbulent jets (Kaminski et al., 2005).In comparison, the second coefficient is still relatively poorly constrained and requires more targeted experiments, both in the laboratory and using three-dimensional numerical models.To capture behaviour in the jet portion of the plume, Kaminski et al. (2005) and Carazzo et al. (2008) developed a modified Reynoldsnumber-dependent entrainment law to account for the negative buoyancy, such that less air is entrained where the plume is denser than the ambient atmosphere.

Modelling plumes in three dimensions
While modifications to one-dimensional axisymmetric models have been relatively minor over the past 60 years, there have been great advancements in the application of more complex three-dimensional models.Great improvements in computational efficiency have enabled the development of increasingly sophisticated models (e.g.Fig. 14d; Cerminara et al., 2016).These models are able to account for a larger range of particle sizes, over much greater scales than possible previously (e.g.Woods, 1988), and are increasingly used to investigate the assumptions utilised in one-dimensional models (Suzuki and Koyaguchi, 2015).Three-dimensional models are increasingly able to account for small-scale processes, for example turbulence and microphysics (Cerminara et al., 2016;Herzog and Graf, 2010;Suzuki and Koyaguchi, 2012), on plume behaviour.Within such three-dimensional models, gas and ash phases are treated as intermingled continua, accounting for mass and momentum transfer between the phases.Despite these great advances in three-dimensional numerical modelling capabilities, there is still a large amount to learn about the relative motion of particles and gas.In addition, further detailed laboratory analysis is required to understand variation in entrainment under real atmospheric conditions, particularly when under the influence of wind.
From a practical point of view, while results are likely more accurate than those from one-dimensional models, the application of three-dimensional numerical models is still limited by computational efficiency, and their results can take many weeks to process and interpret.
12 Ash dispersal Perhaps the most disruptive aspect of an explosive eruption in terms of geographical scale is the injection of volcanic ash into the atmosphere.Ash can be transported hundreds to thousands of kilometres downwind from the source, impacting aviation and downwind communities and infrastructure.Ash dispersal is controlled by a complex relationship between volcanic source and atmospheric conditions; proximity to source dispersion is almost completely controlled by the characteristics of the eruption, while distally, atmospheric physics take over.Given the potential for significant disruption over long timescales, research has focussed on the mechanisms controlling both ash dispersal and deposition.There are a number of similarities between a dispersing plume and the dilute portion of a pyroclastic density current, and therefore many of the scaling laws are the same.
12.1 Laboratory modelling of ash transport processes Laboratory modelling of ash transport in the atmosphere predominantly relates to near-source processes and to investigating controls on sedimentation.Several studies have modelled the volcanic plume as the intrusion of a buoyancy-driven gravity current (e.g.Didden and Maxworthy, 1982;Ivey and Blake, 1985;Bursik et al., 1992;Kotsovinos, 2000) into the atmosphere.In these examples, a dense fluid is injected into a stratified fluid and the less-dense fluid rises through the dense fluid until it reaches its level of neutral buoyancy and begins to spread laterally away from the source.Many experiments have focussed on the physical controls on particle sedimentation within a spreading plume, and constraining parameters such as particle terminal settling velocity.Koyaguchi et al. (2009) considered the effects of turbulence on particle dispersion in the atmosphere by mixing spherical glass-bead particles in water with various intensities of turbulence and measuring the spatial distribution and temporal evolution of the particle concentration.These experiments provide insight into the settling behaviour of particles within a turbulent regime, with results providing information on how particles are dispersed during an eruption.Particle terminal settling velocity of individual particles is estimated by dropping particles with well-constrained characteristics through a fluid of a known density and viscosity (e.g.Dioguardi et al., 2016), with results showing that the size, density, and shape all play a role.
Observations or eruptive events have shown that particle and settling behaviour is often controlled by the characteristics of the mixture rather than individual particles.Laboratory experiments have been used to quantify the effect of the interactions between particles and particle-fluid interaction on sedimentation (Del Bello et al., 2017), by releasing different volume fractions of ash at variable discharge rates through a chamber.Particle settling behaviour was captured using high-speed cameras, which showed a large increase in settling rate with increase in particle volume fraction.The results were validated by a numerical simulation of particle behaviour.Manzella et al. (2015) also looked at the effect of volume fraction of ash on settling behaviour, reproducing the gravitational instabilities noted in field observations by mixing high concentrations of ash into water.An emerging field of interest is exploring how different particles interact with each other, in particular in relation to the formation of particle aggregates (Mueller et al., 2016) where particles are released into a tank and their sticking efficiency is measured using high-speed camera imaging.

Numerically modelling ash dispersion
Dispersion models are used to simulate the dispersal of particulate matter and gases during a volcanic eruption and, in comparison to laboratory models, generally focus on more distal processes.Dispersion models of particulate matter have two main roles: to forecast the dispersion of ash during a volcanic eruption and to reproduce ancient eruptions by fitting model results to observed deposit distributions.Numerical modelling of volcanic ash in the atmosphere requires the definition of three components (Folch, 2012): (1) the source describing the emission of particles and gas in the atmosphere (as explained in Sect.11), (2) an atmospheric model describing the physical characteristics of the environment into which the plume is injected, and (3) the transport model which describes how the particles are transported.
Two main types of numerical models exist: (a) buoyancy models that consider near-source plume characteristics and (b) advection-diffusion ash transport models (Folch, 2012).Buoyant plume models describe the horizontal intrusion of volcanic plumes into the atmosphere as a gravity current (Bursik et al., 1992;Baines et al., 2008;Suzuki and Koyaguchi, 2009;Johnson et al., 2015).They are capable of reproducing both upwind dispersal (e.g.Baines et al., 2008) and plume thickness variation (Johnson et al., 2015), which are key considerations when assessing hazard to aviation.However, the majority of ash models utilised to predict ash transport in the atmosphere are based on the advection of particles by atmospheric winds and the diffusion of ash by atmospheric turbulence.A number of different types of such tephra transport and dispersal models exist, and a comprehensive review is provided in Bonadonna et al. (2012) andFolch (2012).Advection-diffusion models are favoured as they are computationally efficient, allowing results of the order of tens of minutes.Such transport models are typically Eulerian (e.g.TEPHRA2, Bonadonna et al., 2005;FALL3D, Folch et al., 2009;and VAFTAD, Heffter and Stunder, 1993), solving for variable particles at fixed locations, or Lagrangian (e.g.NAME, Jones et al., 2007), calculating the trajectories of a parcel of particles and computing mass concentration by averaging over the background.Hybrid models such as VOL-CALPUFF also exist, which is Eulerian for the plume and Lagrangian in the ash dispersal.Advection-diffusion models predict ash dispersion through the action of vertical wind shear, which can disperse ash in different directions at different altitudes.Their results are therefore highly dependent on the wind field that is used, with offline atmospheric models used in many cases, providing weather information at fixed locations at fixed intervals.
Characterisation of depositional processes is crucial for interpreting volcanic deposits; however these are only ac-counted for in numerical models to a limited extent.Deposition is calculated according to an individual particle's terminal velocity; however observations and laboratory analysis show that the controls on deposition are much more complex, including the formation of aggregates and gravitational instabilities.Recent advances in numerical modelling have aimed to incorporate the effects of aggregation (Folch et al., 2016) and buoyancy forces (Costa et al., 2013) on plume dispersal.Attempts to include aggregation in dispersal simulations are somewhat simplistic and are generally conducted by fixing a priori input grain size distributions rather than accounting for the physics of the process (e.g.Cornell et al., 1983).

Perspectives and conclusions
This review provides an overview of the development of modelling in volcanology, describing some of the first experiments carried out using analogue materials and the development of numerical models to describe volcanic phenomenon.Both laboratory and numerical models are based on theoretical frameworks that have been developed to account for observations and measurements made in nature.These methods have specific challenges that need to be considered when using the modelling results to understand the natural phenomena.Some of the models described above are directly informed by case studies, others are more generalised or focus on a specific process, and where possible the models are tested by comparing the model outputs with expected outcomes based on observed or measured phenomena.In engineering, laboratory models are sometimes used to test and inform the development of numerical models; a numerical model of the laboratory experiment is first created, and the results are then compared with those of the laboratory experiment.This approach has great potential in volcanology, with successful examples evident in studies of lava flow modelling, tephra sedimentation, and plume dynamics; however the integration of laboratory and numerical modelling is yet to be fully explored in volcanology.By assessing the mismatch between laboratory and numerical models and testing the model outputs against natural observations, model errors can be quantified, assumptions investigated, and any limitations in the laboratory or numerical model parameterisation identified.
Our review is relatively unusual in volcanology by considering the insights that have come from both numerical models and scaled laboratory experiments.The scope of our review, which spans from the magma chambers in the lower crust to ash dispersal in the stratosphere, is also vast.However, both of these factors have enabled us to identify emerging lines of thought and make three suggestions below that point towards the future of modelling in volcanology.
Recommendation 1: Utilise a multidisciplinary approach Volcanology as a discipline has significantly benefitted from the understanding and technological advancements across diverse fields, from engineering to material science.Fluid dynamics theory has been applied to study all parts of the volcanic system, for example plume theory (initially developed to study factory emissions) has been applied to the development of ash plumes but also the injection of hot, buoyant magma into the base of a magma chamber.The combination of field observations with geophysical analysis and monitoring techniques has enabled construction of informed conceptual models and hypotheses; these have then been tested in the laboratory using laboratory experiments or computationally using numerical modelling.The limitations of individual techniques are often overcome or mitigated by other methodologies (e.g.Burchardt and Galland, 2016).It is clear that the most significant advancements have come from utilising a multidisciplinary approach, and this needs to be developed further in order to push the frontiers in volcanology.
Recommendation 2: Periodic review of conceptual modelling framework Laboratory and numerical models in volcanology study processes that span the crust-mantle interface into the stratosphere and cross orders of magnitude in time and space.The objectives of any model are carefully defined, and different approaches need to be considered depending on the application.Models that are developed in the laboratory or numerically are ultimately limited by the conceptual models they simulate, and these will be based upon diverse data sources such as geochemistry, petrology, geophysics, real-time observations, and field geology.Referring the model outputs back to known outcomes based on observations from nature is a fundamental step that ensures the model boundary conditions are well informed.
Conceptual models of physical processes in volcanology are constantly evolving, and new insights and developments that come from advancements in formative fields have the potential to revolutionise the framework upon which the laboratory and numerical models are based.For example, discussions on the origin of magma chambers (see Sparks and Cashman, 2017, for a review) and relationship with plutonic bodies (see Lundstrom and Galzner, 2016b, and papers in the same volume for discussions) are being revolutionised by advancements in geochronology as high-resolution dating of ancient magma bodies is now possible, demonstrating that some magma bodies were emplaced over a timescale that is longer than their thermal lifetime (Glazner et al., 2004).Consequently, plutons are now thought to have been accreted by the in situ amalgamation of many small increments and so many numerical and laboratory models that study the dynam-ics of magma chambers and their relationship with plutons need to be revisited.
Recommendation 3: Implement benchmarking, model-intercomparison exercises, and systematic review Model-intercomparison exercises have been employed for modelling tectonic processes (e.g.Schreurs et al., 2006) and similar exercises would be of benefit in volcanology.A recent model-intercomparison of numerical modelling of volcanic plume behaviour was undertaken by Costa et al. (2016) and demonstrated the benefits of this but was limited as several of the models tested have similar underlying assumptions.However, more model-intercomparison exercises are needed for laboratory modelling in volcanology.Systematic evaluation of models has been demonstrated by several studies of volcanic processes, for example in numerical modelling (e.g.Scollo et al., 2008;Costa et al., 2016, and references therein) in which the effect model input uncertainties have on the model outputs, and the interaction of inputs within the model, is evaluated and thus identify important interdependencies.A key advancement in the application of numerical techniques is ensemble modelling, where several different models are applied to investigate the most likely outcomes.Such techniques are well used in other sciences, for example in climate studies, but have yet to be systematically utilised in volcanology.
There is great potential to improve the interaction between laboratory and numerical modelling communities and to develop iterative processes using both techniques to test against data derived from nature.Scaled laboratory experiments allow us to observe and understand physical processes and test basic assumptions; and in parallel with this, numerical simulations model and enable the quantification of processes that are too complex, too large, or last too long to be reproduced in the lab.The approaches are therefore highly complementary, as shown by recent efforts to compare numerical model results and lava flow experiments (Cordonnier et al., 2016 andDietterich et al., 2017) to test the validity of the modelling assumptions.Adopting an engineering approach of systematically testing a numerical simulation against the expected outcomes from a scaled laboratory experiment would be a positive step towards integrating numerical and analogue modelling techniques in volcanology.
We hope this review inspires future, more specialist reviews, for example on vent and conduit processes, which are needed to continue the discussion of laboratory and numerical modelling in different volcanological contexts.

Figure 2 .
Figure 2. Reproduction of Plate IV from Hall (1815).A series of images depict a set of experiments which have since inspired the use of analogue materials to provide a physical explanation for geological observations made in the field.Hall produced an "ideal" coastal section (a) to demonstrate the continuous nature of the folded rock layers observed in the field.Sketches depict a set of experiments that were performed (b) by compressing clay layers to produce convolutions (c) that are reminiscent of the structures observed in the field.

Figure 3 .
Figure 3. Photographs of early analogue experiments that have inspired and informed decades of laboratory studies of volcanic and magmatic processes: (a) excavated plaster of Paris mixture that was injected into a layered tank of gelatine to model dyke and sill emplacement(Hubbert and Willis, 1957), (b) water injected into a freestanding triangular prism of gelatine to simulate magma intrusion in a volcanic rift(Fiske and Jackson, 1972), and (c) injection of lowdensity fluid into a uniform ambient fluid(Morton et al., 1956).

Figure 4 .
Figure4.A flow diagram to represent the optimal approach to use analogue and numerical modelling in volcanology.Observations and measurements from natural volcanic and magmatic phenomena provide the parameters and initial conditions for analogue and numerical models; these models are then tested against nature, with analogue models also aiding the development of numerical models.

Figure 5 .
Figure 5. Flow curves of different fluid rheologies, depending on shear stress and strain rate.

Figure 6 .
Figure 6.Magma rheology studied using analogue materials in laboratory experiments.(a-c) Spherical glass beads, oblate art glitter, and prolate glass fibres in silicone oil (scale bars 1 mm; Mueller et al., 2011); (d) three-phase fluid in which bubbles (black spheres) and spherical glass beads (light translucent particles) are suspended in golden syrup (scale bar 500 µm;Truby et al., 2015).(e) Bubble injection experiments using a small-gap parallel-plate geometry to study the development of permeable pathways in a particle-rich suspension.Particle image velocimetry has been used to measure particle speed in three experiments with different crystal fraction(Oppenheimer et al., 2015).

Figure 7 .
Figure7.Numerical models to simulate magma chamber accretion and associated deformation of the crust: (a) numerical simulation of the magma chamber at Mt Pelée, Martinique, with accumulation of 5 km diameter sills at 15 × 10 −4 km 3 yr −1 over 62 500 years and solved for the temperature of the crust with depth and melt fraction(Annen et al., 2008).(b) Thermo-mechanical model of the magma reservoir at Aira caldera, Japan, using a finite-element model (COMSOL Multiphysics) that incorporates the temperature-dependent viscoelastic rheology of the crust(Hickey et al., 2016).

Figure 8 .
Figure 8. Example numerical models of sheet intrusions: (a) dyke trajectories in the presence of a stress barrier (Maccaferri et al., 2014), (b) the growth of a laccolith(Bunger and Cruden, 2011), and (c) growth of a laccolith or sill with inelastic deformation(Scheibert et al., 2017).

Figure 9 .
Figure 9.A series of photographs demonstrating a range of analogue experiments used to study dyke propagation dynamics.(a) Injection of red-dyed isothermal heptane into a gelatine solid, with l h indicating the length of the dyke's buoyant head(Taisne and Tait, 2011).(b) Injection of a solidifying liquid (wax) causing the formation of an irregular and lobed dyke morphology(Taisne and Tait, 2011).(c) Series of photographs showing the impact of a volcanic edifice on dyke propagation: injection of dyed water into gelatine with a conical surface load of sand and plaster(Kervyn et al., 2009).(d) Schematic sketch and detailed photographs of excavated Vegetaline intrusions such as dykes, hybrids, and cone sheets formed within compacted silica flour experiments(Galland et al., 2014).

Figure 10 .
Figure 10.Analogue experiments studying sill formation.(a-b) Dyke-sill hybrid formation from a feeder dyke in layered gelatine along a weak interface (Kavanagh et al., 2017): (a) polarised light enables stress to be visualised through coloured fringes and (b) a laser-illuminated vertical section through the experiment, showing the sharp boundary between intrusion and host.(c) Digital image correlation of sill formation from a feeder dyke in layered gelatine(Kavanagh et al., 2015).Colours indicate incremental strain in the gelatine, and arrows are displacement vectors at the moment of sill formation: (i) horizontal incremental strain and (ii) vertical incremental strain.(d) Experimental set-up in which Vegetaline is injected into compacted silica flour to model sheet intrusions(Galland et al., 2014).(e-f) Experimental sill formed from a feeder dyke in layered gelatine for which solidification effects are considered(Chanceaux, 2013); (e) typical sill-forming experiment and (f) excavated 3-D morphology of Vegetaline sill showing lobed and segmented propagation front (see alsoChanceaux and Menand, 2014,  2016).

Figure 11 .
Figure 11.Studies of lava dome morphology: (a) analogue experiment by Griffiths and Fink (1997) of polyethylene glycol extruded from a point source imaged in (i) plan view and (ii) side view.(b) Photograph showing (i) external morphology and (ii) cross section through plaster of Paris analogue model of lava dome emplacement seeded with magnetic particles for AMS (scale bars are 10 cm long, photos courtesy of Prokop Zavada.See also Zavada et al., 2009).

Figure 12 .
Figure 12.Analogue models and numerical simulations of lava flows.(a) Effusion of molten wax onto bed of solid wax to study thermal erosion of lava flow into underlying material after (i) 4 and (ii) 14 min (modified from Kerr, 2001).(b) Dehydration of corn starch-water slurry to study the formation of columnar jointing structures in lava flows (modified from Goehring et al., 2006).(c) Cellular automata model of lava flow inundation (green and dark blue) compared with the flow path of the natural lava flow (light blue and dark blue) effused at Mt Etna, Italy (modified from Rongo et al., 2016).

Figure 13 .
Figure 13.(a) Example of early analogue experiments of particulate flows whereby a mixture of methanol, ethylene glycol, and water was released on a slope into a tank of freshwater.The higher density of the mixture in comparison with the ambient fluid means it flows down slope as a gravity current, forming turbulent eddies at the top of the flow(Woods and Bursik, 1994).(b) Small-scale experimental model set-up designed to investigate particle interactions in detail(Roche, 2012).(c) Large-scale analogue experiments using the PELE set-up, allowing large-scale processes within pyroclastic density currents to be investigated.The figure shows emplacement of a dilute mixture of particles and air(Lube et al., 2015).(d) Results from application of the multiphase numerical model PDAC to the blast phase of the Mount St Helens blast of 18 May 1980 (Esposti Ongaro et al., 2012).The model results reproduce numerous features of the flows that originated from the blast, including the formation of turbulent eddies.

Figure 14 .
Figure 14.Examples of numerical simulations and analogue models of volcanic plumes.(a) Flume tank analogue experiments of Carey et al. (1988) in which a particle-laden fluid is injected into another fluid to investigate plume rise and particle sedimentation dynamics.(b) Large-scale experiments to investigate effects of entrainment coefficient on plume rise (Dellino et al., 2014).(c) Modelled representation of a vertical plume using the one-dimensional model PLUME-MoM (de' Michieli Vitturi et al., 2016).(d) Threedimensional model results from the multiphase model ASHEE for the same input conditions as (c) (Cerminara et al., 2016).

Table 1 .
Properties of analogue materials used to model magmas and lavas for rheology and processes.

Table 2 .
Examples of different analogue materials and their combinations used to model different parts of the volcanic and magmatic system.