The use of a shear-thinning polymer as a bubbly magma analogue for scaled laboratory experiments

Analogue materials are commonly used in volcanology to perform scaled laboratory experiments. Analogue ex- periments inform on fundamental ﬂ uid dynamic, structural and mechanical processesthatare typicallyverydif-ﬁ cult to observe and quantify directly in the natural volcanic system. Here we investigate the suitability of an aqueous solution of hydroxyethyl cellulose polymer (HEC) for use as a lava/magma analogue, with a particular focus on its rheological behaviour. We characterize a range of physical properties as functions of the concentra- tion and temperature of the solution: density; speci ﬁ c heat capacity; thermal diffusivity; thermal conductivity; surface tension; as well as rheology. HEC has a non-Newtonian, shear-thinning rheology that depends on the concentration and temperature of the solution. We show that the rheology is well described by the Cross model, which was originally developed for polymer solutions, but has also been applied to bubbly magmas. Using this similarity, an approach for scaling analogue experiments that use shear-thinning polymers, like HEC, to bubbly magma is presented. A detailed work ﬂ ow and a spreadsheet are provided to allow experimentalists to investigate the effects of non-Newtonian behaviour in their existing laboratory set-ups. This contribution willallowforthemorecomplex,butoftenmorerealisticcaseofbubble-bearingmagmastoberigorouslystudied in experimental volcanology.


Introduction
Laboratory experiments offer a way to systematically observe and quantify volcanologically relevant processes that are not accessible to direct observation in the field. Experiments may use natural or analogue materials, and both approaches have advantages and disadvantages. Natural materials, or synthetic materials with the same composition, have the advantage that their behaviour under given physical conditions is expected to replicate directly their behaviour under the same conditions in the volcanic system. However, the experiments must be performed at high temperatures and pressures and, consequently, are usually limited to small sample sizes (Gonnermann et al., 2017;Jones et al., 2016;Quane and Russell, 2005;Renggli et al., 2016;Song et al., 2016). A corollary of the small sample size is that measurements on natural materials tend to be restricted to determining material parameters (e.g. solubility, thermal properties, density, melt viscosity etc.) rather than flow dynamics (Giordano and Dingwell, 2003;Holtz et al., 1995;Lange et al., 1994;Ochs and Lange, 1999;Scheu et al., 2006). Moreover, measurements on natural materials are also generally subject to a high degree of noise due to the natural variability within the samples themselves. This can be somewhat reduced by using synthetic magmas (e.g., Lejeune and Richet, 1995;Yoder Jr and Tilley, 1962).
By contrast, analogue materials are usually chosen to allow observations to take place at temperatures and pressures close to ambient. The focus of many analogue experiments is to study bulk flow dynamics and/or material properties associated with physical interactions among phases. The choice of analogue material is determined by a number of considerations. 1) A specific physical property may be required to facilitate measurements (e.g. transparency for optical measurements). 2) Scaling considerations inform the material properties required for the flow in the natural and experimental systems to be in the same dynamic regime (e.g. laminar vs turbulent, ductile vs brittle, viscousdominated vs surface-tension dominated etc.). 3) It may be desirable to reduce complexity with respect to the natural system in order to isolate and characterize a single process, or a small number of processes, for detailed investigation.
By using scaled, analogue experiments to isolate behaviours and remove multi-component complexity, many fundamental insights have been made to several areas of volcanology. Examples focusing on magma or lava dynamics include: lava flow emplacement (Blake and Bruno, 2000;Castruccio et al., 2010); the rheology of crystal and bubble bearing magmas (Cimarelli et al., 2011;Mueller et al., 2009;Truby et al., 2015); basaltic magma degassing (Namiki and Manga, 2008;Palma et al., 2011) and basaltic eruption dynamics (Del Bello et al., 2015;Jones et al., 2019;Seyfried and Freundt, 2000). Although previous studies have advanced our understanding of volcanic processes, they commonly use Newtonian fluids as analogues. Whilst this simplifies the physics of the problem, in nearly all natural contexts the magma will contain bubbles and therefore exhibit a shear-thinning rheology (Llewellin and Manga, 2005;Mader et al., 2013;Rust and Manga, 2002a). Therefore, to better replicate the natural system, it would be useful to perform experiments using a fluid that has a non-Newtonian rheology under the shear-rates attained. The simplest approach is to add bubbles to a Newtonian analogue fluid (Bagdassarov and Pinkerton, 2004;Rust and Manga, 2002b;Soule and Cashman, 2005). However, the addition of bubbles poses several experimental challenges. Firstly, bubble suspensions are opaque, even at low gas volume fractions. This inhibits visualization of the fluid dynamical processes occurring within the experiment. Secondly, the rheology of the suspension depends on the size distribution of the bubbles, which is difficult to control and measure. Lastly, as bubbles are less dense than the host fluid, they will rise during the experiment adding further complexity to the fluid dynamics operating. We therefore seek an analogue that is transparent and has a shear-thinning rheology that replicates that of a bubble suspension.
Some previous studies have used polymeric solutions for their shear-thinning rheology. To date such studies have been mainly focused around magma degassing and bubble dynamics. Examples include: hydroxypropyl methyl-cellulose to investigate the rise and formation of bubble chains (Kliakhandler, 2002); aqueous solutions of cetyltrimethylammonium bromide and sodium salicylate to study degassing dynamics through a thin viscoelastic layer ; a commercial hair dressing gel to investigate magma degassing regimes (Divoux et al., , 2009; aqueous solutions of polyacrylamide and carboxymethylcellulose (CMC) mixed with glycerol and water to investigate stress and relaxation around bubbles during their rise (Li et al., 2004); aqueous solutions of CMC in exchange flow experiments (Huppert and Hallworth, 2007); aqueous solutions of celluloseester in analogue experiments for basaltic conduit flow (Seyfried and Freundt, 2000); and aqueous solutions of hydroxyethyl cellulose (same chemical name as used in this study but sourced from a different manufacturer with an unknown molecular weight) in bubble growth and expansion experiments (Zhang, 1998;Zhang et al., 1997). Although some of these previous studies describe the rheology of these non-Newtonian fluids as a function of strain-rate, scaling to the natural system has been limited to simple comparison of rheological properties (e.g. both magma and the chosen analogue fluid have a shear-thinning rheology; Divoux et al., 2011Divoux et al., , 2009Vidal et al., 2011). Consequently, there is no complete and rigorous scaling to the natural volcanic system.
Here, we characterize a shear-thinning polymer (hydroxyethyl cellulose, HEC) that is viscoelastic, transparent and has a highly adjustable viscosity over a commonly used shear-rate range making it a suitable analogue material for many existing scaled experimental set-ups. Specifically, we use HEC manufactured by the Dow Chemical Company with the trade name Cellosize™ and grade QP 52000H. We focus on quantifying the polymer rheology, since that is essential for scaling any analogue experiment in fluid dynamics. We also determine the density and thermophysical properties (specific heat capacity, thermal diffusivity and thermal conductivity), which are important for scaling experiments involving convection, and surface tension, for experiments involving bubble or gas slug rise, for example. The discussion focusses on the use of HEC and shear-thinning polymers in scaled experiments of volcanic flows. We find that HEC provides a useful analogue for bubbly magma over the conditions commonly encountered in existing experimental devices. For a full list of parameters and associated notation used throughout this paper, refer to Table 1.

Chemical properties
Hydroxyethyl cellulose (HEC) is a water-soluble polymer that is used in the cosmetics, food, and pharmaceutical industries as a thickening and gelling agent. The chemical structure is shown by Fig. S1 in the Supplementary Information. HEC is available under various trade names, including Natrosol™ (Ashland Inc.), Clearcel™ (Chemcolloids Ltd.), and Cellosize™ (Dow Chemical Company). The latter is used in this study. Cellosize™ is available in two types: QP which disperses rapidly and dissolves in solution, and WP which hydrates rapidly and tends to aggregate. Each type comes in a range of viscosity grades (Dow, 2005). In this study, we use grade QP 52000H because it dissolves easily and, at low shear rates, has a viscosity similar to commonly used analogue fluids (e.g. golden syrup). In general, higher molecular weight grades of HEC have higher viscosity at the same concentration, and are more strongly shear thinning. This behaviour arises principally from the greater degree of entanglement of higher molecular weight polymer chains.

Preparing HEC solutions
Solutions of 0.5, 0.75, 1, 1.25, 1.5 and 1.75 wt% concentrations (X c ) of Cellosize QP 5200H (hereafter referred to as 'HEC') in tap water (pH 7.85 ± 0.03) were produced. For a 1 wt% solution we add 1 g of HEC powder to 99 g of water. HEC is hygroscopic and therefore gains up to 2.5% mass by the addition of moisture at standard laboratory conditions (5-6 g m −3 humidity); in this study we use HEC that has been stored and prepared under these typical conditions. The uncertain degree of hydration leads to an uncertainty in the weighed mass of the HEC powder of up to 2.5%; hence, the uncertainty in the concentration of the resulting solutions is also approximately 2.5% of the calculated valuethis uncertainty should be propagated when applying any of the empirical models for material properties presented later in Sections 3 and 4. The uncertainty resulting from variable hydration dominates the small error associated with the precision of the balance used to weigh the powder and water. The beaker containing the solution was placed in a bowl of warm water (~55°C) and stirred slowly, to avoid the introduction of bubbles, for 2-10 min until it thickened. The time for thickening depends on many factors such as temperature (higher temperatures reduce thickening time) and pH (higher pH values reduce thickening time). The thickened solution was transferred into an airtight container for storage to minimize dehydration, contamination and mould formation. The rheology of HEC solutions can depend on the pH of the water, so it is important, when preparing HEC solutions, to ensure that a stable pH is maintained (Dow, 2005).
All the samples characterized in this study were of small volume (~300 ml). If large volumes (of several litres or more) of HEC are required, then a different approach is recommended. Firstly, put a known volume of hot (~50°C) water in a large mixing bucket and slowly and steadily, over the course of a few minutes, add the HEC powder to reach the desired concentration. During the gradual powder addition, the solution should be stirred continuously with an electric paddle mixer to avoid aggregation until the solution thickens (around 10 min, but longer for colder water). Given the greater uncertainty in the concentration value when preparing large volumes we recommend that the concentration should be determined by characterising the rheology of an aliquot.

Magma rheology
Pure silicate melts are Newtonian for strain-rates Á γ≪1=λ r where λ r [s] is the structural relaxation time of the melt Webb, 1990, 1989;Giordano et al., 2008;Giordano and Dingwell, 2003;Hess and Dingwell, 1996). A Newtonian liquid deforms with strain-rate γ • [s −1 ] proportional to the applied stress τ [Pa]; the constant of proportionality is the viscosity μ 0 [Pa s]: The viscosity of natural silicate melts spans a huge range,~10 −1 -10 12 Pa s, as silica content, temperature and dissolved water content all exert a strong control on melt viscosity (Giordano et al., 2008). The viscosity of silicate melts is comprehensively reviewed by Giordano et al. (2008).
In volcanic systems pure melts are rare, and magma is usually a multiphase suspension composed of a melt (the liquid phase with Newtonian viscosity μ 0 ) that suspends bubbles (gas) and/or crystals (solid).
The addition of bubbles and/or crystals induces a strain-rate dependent (non-Newtonian) rheology. It is therefore useful to define an apparent viscosity η [Pa s], which is the ratio of stress to strain-rate at any given strain-rate. Furthermore, the ratio of the apparent viscosity of the suspension to the viscosity of the suspending liquid is termed the relative viscosity η r = η/μ 0 . In this work we focus on the rheology of bubble suspensions. Bubbles introduce shear-thinning rheology that can be characterized by the dimensionless capillary number Ca.
where λ b = μ 0 a/Γ is the bubble relaxation time [s], a [m] is the undeformed spherical bubble radius, and Γ [N m −1 ] is the surface tension. At low Ca the surface tension stress dominates so the bubbles remain close to spherical in shape and act to increase the suspension viscosity. Conversely, at high Ca viscous stress dominates so the bubbles become deformed and act to decrease the suspension viscosity (Llewellin and Manga, 2005;Rust and Manga, 2002a;Spera, 2002, 1992).
Previous rheological measurements on both bubbly analogue suspensions and bubbly silicate melts have shown that the rheology of a bubble suspension is shear thinning (Llewellin et al., 2002;Rust and Manga, 2002a;Stein and Spera, 2002;Truby et al., 2015). To illustrate this we have compiled data from three of these previous studies that use suspensions with several different bubble volume fractions ϕ deformed at a range of capillary numbers (Fig. 1a). These data show that the addition of a suspended gas phase to a Newtonian liquid introduces non-Newtonian behaviour, such that Eq. (1) is insufficient to describe the rheology. At high and low Ca, asymptotic limits of relative viscosity are reached (Llewellin and Manga, 2005;Pal, 2003). Mader et al. (2013) present these limits of zero-shear relative viscosity η r0 and infiniteshear relative viscosity η r∞ as functions of bubble volume fraction: where η o and η ∞ are the apparent viscosity of the suspension in, respectively, the low and high strain-rate limits: At intermediate shear rates, the relative viscosity falls between η r∞ and η r0 , and is a function of Ca as well as ϕ. The viscosity can be written in the form of a Cross model valid for steady simple-shearing flow (Mader et al., 2013;Rust and Manga, 2002b): where K = 6/5 and m = 2 for a monodisperse bubble population. Polydisperse suspensions can be accounted for by adjusting m to lower values of~1 (Mader et al., 2013). The Cross model for bubble suspensions assumes laminar flow and has only been validated on suspensions with gas volume fractions in the range 0 ≤ ϕ ≤ 0.46. Upon transition from a bubble suspension to a polyhedral foam, which occurs at ϕ~0.74, the rheological characteristics are likely to be very different (Mader et al., 2013). Our analysis presented herein therefore assumes laminar flow and a bubble suspension, rather than foam, rheology.

Polymer rheology
Some polymers also show shear-thinning rheology (Fig. 1b). As with bubble suspensions, these polymers have an apparent viscosity that tends to an asymptotic value at high and low values of strain-rate, as given by Eqs. (5) and (6). The Cross model (Cross, 1965) was originally developed to describe polymer rheology and is widely used in polymer science, where it is written in the form (compare with Eq. (7) above): where η is the apparent viscosity of the polymer [Pa s], c [s] is the Cross constant and p is the Cross exponent (Fig. 1b). Increasing the Cross exponent increases the gradient of the curve between the viscosity limits, whilst varying the Cross constant translates curves along the strain-rate axis (Fig. 1b). When using Eqs. (3)-(8) particular care should be taken with subscripts to ensure that the correct quantity is used. The subscripts denote different viscosity values (e.g. apparent vs. relative viscosity) and, for the reader's ease, a list of notation is provided in Table 1.
The fact that the rheology of both HEC and bubble suspensions can be described using a Cross model makes it straightforward to scale the rheologies of the two fluids. This is explored in detail later, in Section 5.

Viscoelasticity and the Cox-Merz rule
Both bubble suspensions and polymer solutions can be viscoelastic. This means that the deformation in response to a driving stress can involve both elastic (recoverable) and viscous (non-recoverable) strain. The degree of viscous and elastic response depends on the timescale of deformation in relation to the relaxation timescale of the material (Barnes et al., 1989). Viscoelastic behaviour can be investigated under oscillatory rheometry (Section 3.4) in which a sinusoidally varying stress is applied to the sample.
The complex shear modulus G * [Pa] describes the resistance to deformation of a viscoelastic material. G * is comprised of the storage modulus G′ [Pa] and the loss modulus G″ [Pa] which characterize elastic and viscous behaviour respectively: . Under oscillatory rheometry, frequencies can be identified for which viscous (G″ N G′) or elastic (G′ N G″) behaviour dominates (Ferry, 1980). The phase shift δ between the applied stress waveform, and resulting strain waveform describes the relative contribution of viscous and elastic deformation such that δ = 0°for purely elastic materials and δ = 90°for purely viscous (Newtonian) materials. Furthermore, the Maxwell relaxation time λ m [s] can be defined as the inverse of the frequency at which G″ = G′ in the frequency [Hz] domain (Schellart, 2011).
It is also common for viscoelastic materials to be described by the complex viscosity η * [Pa s], rather than the complex shear modulus: where η′ [Pa s] is the dynamic viscosity and η″ [Pa s] is the storage viscosity. These parameters are related to the complex shear modulus via the angular frequency of oscillation ω [rad s −1 ] such that η * = G * /iω. A purely empirical relationship described by Cox and Merz (1958) has provided a link between steady state and oscillatory rheometry measurements. They observed ηðγ • Þ ¼ jη Ã ðωÞj, i.e. that the magnitude of the complex viscosity as a function of angular frequency under oscillatory rheometry is the same as the apparent viscosity as a function of strain-rate under steady rheometry; this defines the Cox-Merz rule (Cox and Merz, 1958). The Cox-Merz rule has been shown to be approximately valid for bubble suspensions (Mader et al., 2013). Flow conditions are unsteady under oscillatory rheometry, so the capillary number Ca is undefined, hence Eq. (7) cannot be applied. Llewellin et al. (2002) consider flow of bubble suspensions under unsteady conditions, and introduce the dynamic capillary number is the timederivative of the strain-rate. Under oscillatory rheometry, this becomes Cd = λ b ω, and the Cox-Merz rule can be expressed as: This approximation underpins a more-general form of Eq. (7) that can be applied to both steady and unsteady bubbly flows, under arbitrarily large strains (Mader et al., 2013):

Method of measurement
Two types of rheometry were performed on HEC solutions; firstly, rotational rheometry was used to characterize the steady-state viscosity and shear-thinning rheology. Secondly, oscillatory rheometry was used to quantify both the viscous and elastic components of deformation (Mezger, 2006;Whorlow, 1992). All the rheological characterization was performed on a TA Instruments Discovery Hybrid HR2 rheometer at Durham University, using a cone and plate measuring geometry, with a temperature-controlled Peltier plate (Fig. 2). In all experiments a solvent trap, filled with deionised water, was used to prevent the solutions dehydrating during measurement. Before the rheometric measurements, the rheometer and sensor system were calibrated for internal friction, thermal expansion, inertia and precession upon rotation. Each time a new aliquot of sample was loaded, it was left for 10 min to reach thermal equilibrium.
For the rotational rheometry the shear stress was measured at strain-rates from 0.01 s −1 to 1000 s −1 with 10 measurement points per decade, distributed logarithmically. Using this experimental approach, viscosities were determined for all HEC concentrations at temperatures: 5°C, 10°C, 15°C, 20°C, 30°C, 40°C, 50°C and 60°C. Lastly, to test HEC stability over time, a 1.25 wt% solution was measured at 20°C after 0, 9, 11, 15 and 21 days stored in an air tight container at lab temperature.
For the oscillatory measurements, a stress amplitude sweep (0.01 to 1000 Pa) was first performed at a 10 rad s −1 constant angular frequency. From this test the material's linear viscoelastic region (LVE) was defined; for stresses in this region material properties such as the complex viscosity, phase angle, loss and storage modulus are approximately constant. A frequency sweep (0.1 to 100 rad s −1 ) was then performed at a constant stress within the LVE region (1 to 0.1 Pa in this Cone and Peltier plate measuring geometry used with the TA instruments Discovery Hybrid HR2 rotational rheometer. After the sample was loaded the solvent trap was filled with deionised water during each experiment to prevent sample dehydration upon heating.  8)). Error bars are smaller than data points; maximum 1σ = 0.77 Pa s calculated with seven repeat measurements when deliberately slightly under-and over-filling the sample. study) in order to characterize the time-dependent rheology. All oscillatory tests were performed at 20°C.

Results and analysis
3.5.1. Shear-rate dependent viscosity All concentrations of HEC exhibit a shear-thinning rheology (Fig. 3) that is well described by the Cross model (Eq. (8)) with η 0 , η ∞ , c, and p as fitting parameters. For shear rates in the measurable range (0.01 s −1 to 1000 s −1 ) the data tend toward a zero-shear apparent viscosity η 0 at low strain-rates, and a much lower infinite-shear apparent viscosity η ∞ at high strain-rates. Measurements become unreliable at the extremely low viscosities encountered at very high strain-rates, so we are not able to reach the infinite-shear apparent viscosity η ∞ . Nonetheless, the Cross model (Eq. (8)) provides a good fit to the data across the measured range. Curve fitting is subject to the constraint that η ∞ must be greater than or equal to the viscosity of water at the experimental temperature. This prevents the production of physically-implausible values of η ∞ during model fitting; in practice, the best-fit value for η ∞ is always equal to the viscosity of water (Table 2). Comparing HEC concentrations ( Fig. 3a to f) it can be seen that solution viscosity increases with increasing concentration. Also, for all concentrations, increasing the temperature reduces the viscosity of the solution at all strain-rates, although this reduction is most pronounced at low strain-rates. A full list of Cross model parameters derived from fitting to data are presented in Table 2.
The systematic changes in Cross model fits shown in Fig. 3 are encapsulated in systematic variation in the parameters η 0 , c and p with temperature and concentration (cf. Fig. 1b). For a given HEC concentration these Cross parameters (η 0 , c and p) can be expressed as functions of temperature in the form: η 0 = A 1 e A2T ; c = B 1 e B2T and p = C 1 T + C 2 , where T is in Celsius. These relationships are purely empirical and in each case we chose the simplest functional form that gives a good fit to data with two fitting parameters. Fitting parameter values can be found in Table 3.
By adding the Cross parameters' dependency on HEC concentration we produce a set of purely empirical Eqs. (13)-(15) for η 0 , c and p as a function of temperature T [°C] and concentration X c . These fits are shown by the solid lines in Fig. 4. All the fits presented here are for HEC concentrations ≥0.75 wt%, the 0.5 wt% data set has been excluded from these fits as values are close to the rheometer sensor system sensitivity limit (Fig. 4) and show considerably more scatter than for the higher concentrations. The temperature and concentration dependence of the zero-shear apparent viscosity (Fig. 4a) is well reproduced using an exponential equation: where D 1 = 0.663 Pa s, D 2 = 3.99 and D 3 = 0.0516°C −1 ; hence, η 0 increases exponentially with HEC concentration and decreases exponentially with temperature. The only previous work describing the rheology of HEC that we have found (Haas and Eloranta, 1965) characterizes the zero-shear apparent viscosity for seven different concentrations at 25°C, and three different concentrations at 48°C. They propose a power-law relationship between zero-shear apparent viscosity and concentration at 25°C. The values of η 0 that they report are roughly two orders of magnitude higher than those in our study, at . Note that η ∞ is simply the viscosity of water at the specified temperature. These individual models for each concentration are more accurate than the more general models (Table 3). the same concentration, suggesting that they used a HEC formulation with a different molecular weight that is unspecified in their paper. The Cross constant c can also be described using an exponential form; c decreases rapidly with increased temperature such that: where E 1 = 11.1 s, E 2 = 7.00 s and E 3 = −0.0306°C −1 . Finally, the Cross exponent p can be best described by a linear relationship: where F 1 = 0.0973, F 2 = 0.00217°C −1 and F 3 = 0.615. We also report changes in the solution apparent viscosity through time after the initial preparation of the solution (Fig. 5). It should be noted that the 500 ml sample (from which measured aliquots were taken) was stored in an airtight container at standard lab temperatures (18-21°C). In general, the tests show a progressive decrease in viscosity with age which is most pronounced at low-shear (Fig. 5b). If the samples were left open to contaminants, additional and more rapid degradation would be expected. Qualitative observations suggest that additional dehydration of HEC in an open container would generate a high viscosity layer on the surface in contact with air. We suggest that HEC solutions should be prepared immediately before use, particularly in experimental set-ups using low shear-rates, and should not be stored for long periods of time.  (15) respectively. To determine the error associated with Cross model fitting 10,000 least square regressions were applied to the raw rheometry data, randomized within their 1σ range. This gave 1σ values of 0.00047, 0.00047 and 0.00002 for η 0 , c and p respectively. Error bars are smaller than the data points. Data scatter is primarily associated with poor instrument sensitivity at low viscosities.

Viscoelasticity
From the amplitude sweep (Fig. 6a) we can identify the LVE region; this region extends to higher stresses with increased HEC concentration. Stresses lower than 1 Pa (for X c ≥ 0.75) or 0.1 Pa (for X c b 0.75) fall within the LVE region. Over the frequency range measured (0.1 to 100 rad s −1 or 0.016 to 15.92 Hz) all the HEC concentrations display viscoelastic behaviour (0°b δ b 90°). Elastic behaviour increases with increasing solution concentration and with increasing oscillation frequency (Fig. 6b). All the HEC solution concentrations except from 0.5 wt% show an intersection between G′ and G″ (Fig. 6c). The position of this intersection in the frequency [Hz] domain gives the Maxwell relaxation times λ m [s] (Table 4). We also present the complex viscosity as a function of frequency for the HEC solutions studied (Fig. 6d). By comparison with the rotational rheometry measurements (lines in Fig. 6d) it is clear that HEC solutions follow the Cox-Merz rule given by Eq. (11) to a close approximation.

Other properties of HEC
In this section, we present data for a range of other physical properties that are relevant to the effective scaling of HEC to magma in the volcanic system. These properties include density, thermal conductivity, thermal diffusivity, specific heat capacity and surface tension. For each physical property we describe the volcanological motivation, the measurement technique, and typical values for the HEC polymer. We note that HEC has photoelastic properties (Dow, 2005) this is not investigated further here, but could potentially be exploited to visualize stresses within the material during flow.

Background and magmatic properties
Measurements on multi-component silicate liquids have shown that their density is linearly dependent on temperature and also depends on oxide concentration, such that the addition of low density oxides (e.g., K 2 O and Al 2 O 3 ) act to lower the silicate liquid density and the addition of high density oxides (e.g., MgO and FeO) act to raise the liquid silicate density (Lange and Carmichael, 1987;Lesher and Spera, 2015). It is also known that small concentrations of dissolved water and carbon dioxide have a large effect on reducing the silicate melt density. For example, adding 1 wt% water to a dry basaltic melt reduces its density bỹ 70 kg m −3 (Lange and Carmichael, 1990). Density is a parameter that often enters into the scaling of fluid dynamic studies; one example is the Reynolds number (Re), defined as: Re = ρνl/η. Where ρ [kg m −3 ] is the density, ν [m s −1 ] is the fluid velocity, l [m] is a characteristic length scale, and η [Pa s] is the apparent viscosity. The Reynolds number is the ratio of inertial to viscous forces and is commonly used in volcanology, and beyond, to distinguish between laminar and turbulent flow regimes. In analogue experiments we often wish to match the Reynolds number of our experiments to the natural system. For a given experimental set-up, the fluid viscosity and density must be adjusted to give an appropriate Re. In most cases it is more practical to change the viscosity and measure the fluid density.

Method of measurement
HEC solutions were added to a calibrated 50 ml pycnometer flask, heated to the measurement temperature and left to equilibrate at 10,  8)) using parameters derived from steady rheometry (Fig. 4, Eqs. (13)-(15)). The close fit to the oscillatory data confirms that the Cox-Merz rule holds for these solutions. Error bars on all subplots are smaller than the data points. 20, 30, 40, 50, 60 and 70°C within a temperature-controlled water bath. During heating the HEC solution was able to overflow from the pycnometer flask. Once equilibrated, the flask exterior was dried and then its mass measured on a ±500 μg desktop balance.

Results and analysis
The density of HEC solutions is found to vary both as a function of temperature and concentration (Fig. 7). Over the temperature range investigated (10 to 70°C) the density of a particular HEC concentration follows the purely empirical general expression: where T [°C] is the temperature and G 1 [kg m −3°C−2 ], G 2 [kg m −3°C−1 ] and G 3 [kg m −3 ] are constants ( Table 5). The density measurements shown in Fig. 7 have a maximum uncertainty of ±0.168 kg m −3 . Furthermore, the concentration dependency at a specific temperature follows the empirical form: where H 1 [kg m −3 ] and H 2 [kg m −3 ] are constants reported in Table 6. HEC density is very similar to water, we therefore suggest that HEC will be particularly useful for experiments that commonly use water and in which buoyancy is important (e.g. analogue magma intrusion experiments; Kavanagh et al., 2018).

Background and magmatic properties
Thermophysical parameters (specific heat capacity, thermal diffusivity, and thermal conductivity) are important in a wide range of volcanological problems (e.g. magma convection, cooling of lava flows). The specific heat capacity (C P [J K −1 kg −1 ]) is the amount of energy required to change the temperature of a unit mass of a substance by one Kelvin: where Q [J] is the amount of heat added, m l [kg] is the mass of liquid and T [K] is the temperature (in Kelvin for the thermophysical parameters, unlike all previous cases where T was in°C). Specific heat capacity of mafic to intermediate magmas commonly falls within the range of~60 to 1500 J K −1 kg −1 (Bachmann and Bergantz, 2006;Bohrson and Spera, 2001;Di Genova et al., 2014;Dobran, 2012;Robert et al., 2015Robert et al., , 2014Spera, 2000;Spera and Bohrson, 2001). Giordano and Russell (2016) present an empirical model for the specific heat capacity of silicate melt and show that it increases linearly with silica content. Thermal diffusivity α [m 2 s −1 ] is a measure of the rate of conductive heat transfer within a material. A thermal diffusivity of 10 −5 m 2 s −1 has been reported for a basalt (Huppert and Sparks, 1988; Jellinek and Kerr,   Table 5. 1σ error bars are smaller than the data points. (b) The linear relationship between density and concentration, colour contoured for temperature. Fits in the form of Eq. (17) using the values listed in Table 6.

Table 5
Density fitting parameters as a function of temperature (shown in Fig. 7a).  Table 6 Density fitting parameters as a function of concentration (shown in Fig. 7b). 1999), whereas a thermal diffusivity of 10 −6 m 2 s −1 has been considered typical for a more silicic melt (Bachmann and Bergantz, 2006). The compositional dependency of silicate melt thermal diffusivity is comprehensively reported in Hofmeister et al. (2016). Models for thermal diffusivity as a function of silicate melt temperature can be found in Bagdassarov and Dingwell (1994), Romine et al. (2012), Hofmeister et al. (2016 and Wadsworth et al. (2017), for example. Thermal transport may alternatively be encapsulated in the thermal conductivity κ [W m −1 K −1 ], which can be expressed as a combination of parameters presented previously: This parameter has been previously used to numerically model the effects of magma injection on the thermal structure of the crust (Annen and Sparks, 2002). A thermal conductivity of approximately 2-5 W m −1 K −1 is considered typical for a silicate melt (Bachmann and Bergantz, 2006;Bartlett, 1969;Chapman and Furlong, 1992;Dobran, 2012).

Method of measurement
The specific heat capacity C p (Eq. (18)) was measured using a NETZSCH 204 F1 Phoenix differential scanning calorimeter (DSC) by NETZSCH instruments. Values reported for all HEC concentrations are at 25°C. Additionally, the 1.75 wt% sample is reported at 10, 30, 40 and 50°C to investigate the temperature dependence of specific heat capacity. Each test used a~10 mg aliquot contained within an aluminium crucible. All tests were conducted in a dynamic nitrogen atmosphere (20 ml min −1 ) and heated from −5°C to 55°C at a rate of 10 K min −1 .
The thermal diffusivity α of each HEC concentration was measured using a NETZSCH laser flash apparatus (LFA) 467 Hyperflash. During the laser flash technique, the front of the sample was heated by a xenon flash lamp and the resultant temperature increase on a parallel rear plane was measured using an IR-detector. The thermal diffusivity is then calculated by the following relationship: α = 0.1388x 2 /t 1/2 where x [m] is the sample thickness and t 1/2 [s] is the time taken for the sample to reach half of the temperature maximum. Each sample was measured five times and averaged before the LFA moved to the next temperature step. Lastly, the thermal conductivity κ was calculated via Eq. (19) for each HEC concentration, using measured values of ρ, C p and α. The errors for C p and α are derived from the maximum standard deviation reported by NETZSCH instruments and then compounded for the error associated with κ. Pure water values were gathered from the international steam tables (Grigull et al., 2012) for ρ, C p and κ and used to calculate α. These values were not used in the model fits because they do not represent the exact water used for preparing the HEC solutions, but are included as a guide.

Results and analysis
The three thermophysical properties C p , κ and α are presented as a function of HEC concentration (Fig. 8). None of the thermophysical parameters show a dependence on temperature within the measured range 10 to 70°C (Fig. S2 in Supplementary Information). For HEC solutions at 25°C, C p decreases with increasing concentration. Over the concentration range investigated this relationship is linear and can be described by: where I 1 = 77.38 J K −1 kg −1 and I 2 = 4201 J K −1 kg −1 . The thermal diffusivity can also be described by a linear correlation (Fig. 8b) and increases with increasing HEC concentration such that: where J 1 = 1.33 × 10 −9 m 2 s −1 kg −1 and J 2 = 1.42 × 10 −7 m 2 s −1 . Lastly, the thermal conductivity is found to show no systematic variation with HEC concentration (Fig. 8c). The mean κ is 0.592 W m −1 K −1 for the concentration range 0.5 to 1.75 wt%. Both the thermal diffusivity and conductivity of HEC are at least one order of magnitude lower than magmas. This is desirable for scaling of analogue experiments, in which length scales are typically shorter than in the natural system.

Background and magmatic properties
In general, surface tension forces act to reduce surface area, for instance by restoring a bubble or melt droplet to a spherical shape. Surface tension is an important parameter that forms a part of many dimensionless groups that are commonly used in experimental scaling, such as the Morton number Mo ¼ gμ 0 4 =ρΓ 3 , which can be thought of as a material property group that captures the balance of stresses arising from viscosity, inertia, and interfacial tension, and the Eötvös number (or Bond number) Eo = ρgd 2 /Γ, where d [m] is the bubble or droplet diameter, which is the ratio of stresses arising from buoyancy and interfacial tension. These groups are particularly important when performing analogue experiments on bubble/gas-slug rise and magma degassing (e.g., Divoux et al., 2009;Kliakhandler, 2002;Llewellin et al., 2011;Seyfried and Freundt, 2000). In such experiments, HEC could be used as a transparent analogue for magma containing small, coupled bubbles, through which large, decoupled bubbles or slugs are rising, allowing the behaviour of the large bubbles to be observed directly.

Method of measurement
The surface tension measurements were conducted by the Dow Chemical Company on 0.01, 0.1 and 1 wt% HEC solutions at 25°C following the ASTM D1331 standard and reproduced here from the Cellosize™ product data booklet (Dow, 2005). The method involves a calibrated tensiometer with an attached, clean platinum ring positioned above the sample vessel. The sample is slowly raised until the platinum ring is submerged, then slowly lowered whilst tension is applied to the wire to keep the suspended ring in a fixed position. As the ring detaches from the liquid surface, the dial reading on the tensiometer is recorded. The dial reading is then converted into a true surface tension by a correction factor (Harkins and Jordan, 1930).

Results and analysis
Surface tension values reported in Table 7 are slightly lower than for pure water (~72 mN m −1 at 25°C; Vargaftik et al., 1983), and do not show a strong dependence on HEC concentration. A maximum variation of 0.7 mN m −1 is observed (roughly 1% of the mean value) when changing an order of magnitude in solution concentration.

Discussion
Scaling analogue experiments that use Newtonian fluids is straightforward because the flow behaviour is described by a single parameter, the Newtonian viscosity. For non-Newtonian fluids it is less straightforward. In this section, we demonstrate how HEC can be scaled as an analogue for bubble suspensions, which have a shear thinning rheology. In this case, the scaling involves the (Newtonian) melt viscosity, the gas volume fraction, which controls the magnitude of the viscosity change due to shear thinning, and the relaxation time, which controls the strain-rate dependence of viscosity.

Comparing a shear-thinning polymer to a bubble suspension
The rheology of both bubble suspensions and HEC can be described using a Cross model (Eqs. (7) and (8) respectively), re-written below as Eqs. (22) and (23) for the reader's ease: The model for bubble suspensions is in dimensionless form, so, to facilitate comparison of the two models, we put the model for HEC in dimensionless form too. By analogy with the bubble suspension rheology case (Eqs. (3) and (4)), we non-dimensionalize the low and high strainrate asymptotic viscosities by a nominal Newtonian viscosity μ n where η n0 and η n∞ are the dimensionless zero-shear relative viscosity and dimensionless infinite-shear relative viscosity respectively. We also define the dimensionless 'Cross number' N C , such that Eq. (23) becomes: where η n is the dimensionless relative viscosity. Our goal is to identify a HEC solution that has approximately the same dimensionless rheological behaviour as a bubble suspension of interest, such that η n (N C ) ≈ η r (KCa). If such an approximation can be found, then there is equivalence between Eqs. (24) and (3) and Eqs. (25) and (4) η n0 ≡ η r0 ð27Þ and moreover where ≡ denotes an equivalence between the parameters. These parameters have the same dimensions and if perfectly scaled would be equal. Inspection of the shapes of the η n (N C ) curves for HEC and η r (Ca) curves for bubble suspensions (Fig. 9a gives an example) shows that HEC always has a gentler slope linking the high and low viscosity asymptotes than a bubble suspension (i.e. p b m), and the low viscosity asymptote is always lower than for a bubble suspension (i.e. η n∞ b η r∞ ). Consequently, it is not possible to have a perfect rheological analogy between the two materials, and so the parameters in Eqs. (28) and (30) cannot be matched exactly. Nonetheless, it is possible to obtain a good approximate agreement between rheologies over a limited range of Ca and N C , particularly for Ca b 10 in this example. Note that in Fig. 9a the viscosity of the HEC solution is normalized so that Eq. (26) is satisfied. From Eqs. (24) and (27) we have μ n = η 0 /η r0 hence, from Eq. (3), it follows that μ n = η 0 /(1 − ϕ) −1 .
We can improve on the analogy shown in Fig. 9a by restricting our attention to lower capillary/Cross numbers, and scaling the Cross number to shift the HEC curve toward the bubble suspension curve (Fig. 9b).
To do this, we define the scaled Cross number as N Cs ≡ N C /b, and replace N C with N Cs in Eq. (26).
We then adjust the dimensionless 'scale factor' b so that there is a good fit between η n (N Cs ) and η r (KCa) over the Ca range of interest. So an improved fit is produced if we put (cf Eq. (29)): Table 7 Surface tension values for a range of HEC concentrations. Data taken from Dow chemical hydroxyethyl cellulose data sheet (Dow, 2005 It is now possible to arrive at a scaling between the strain-rate in the HEC, and the effective capillary number of the bubble suspension for which it is an analogue: We term the group bc/K the 'strain-rate factor'. Similarly, the viscosity of the HEC solution scales to the viscosity of the bubble suspension via the 'viscosity ratio':

Scaling an analogue experiment: an example calculation
For the users ease a spreadsheet with supporting text that simplifies the scaling for shear thinning rheology is provided in the Supplementary Information. The spreadsheet details how, by adjusting the HEC solution concentration and temperature, it is possible to replicate a bubbly magma with a specified gas fraction and melt viscosity, over a specified Ca range. Now, as an example, we consider scaling analogue experiments using HEC to investigate the flow of bubbly basaltic magma along a dyke of width 1 m. Suppose that the magma has the following properties: the pure (bubble-free) melt has a viscosity of μ 0 = 1000 Pa s; the gas volume fraction ϕ = 0.4; the Cross model exponent m = 1.2 (i.e. the suspension is polydisperse); the mean bubble radius a = 1 mm; and the melt-gas surface tension Γ = 0.2 N m −1 . We wish to investigate the role that shear thinning plays on flow organization, up to Ca = 10. This requires that we can find a HEC solution that scales to the appropriate rheology over the range 0 b Ca ≤ 10, and that falls in the same dynamic regime (i.e. Reynolds number is similar).
First, we estimate the Reynolds number of the natural flow. The maximum capillarity is reached at the dyke wall. To be conservative in our estimation of Reynolds number, we assume that, when Ca = 10 at the dyke wall, the bulk of the magma is in the high capillary number regime, and that we can therefore treat the flow as if the magma were Newtonian (Fig. 1), with apparent viscosity equal to the low viscosity asymptotic value for the bubbly magma (i.e. η = μ 0 η r∞ = 426 Pa s). We further assume that the velocity profile is parabolic (i.e. we are at sufficiently low Reynolds number for Poiseuille flow) in which case the strain-rate at the wall is related to the driving pressure gradient ΔP/L according to where D is the width of the dyke. The volume flux q per unit length of dyke is given by and the mean velocity of the flow is given by ν = q/D. The Reynolds number is given by Re = ρDν/η. Evaluating these Equations for the current case gives Re = 1.3. Our assumption of Poiseuille flow and low Re is therefore valid and our Cross model-based scaling approach can be used. Note that the Newtonian and parabolic assumptions are purely for the purpose of conservative estimation of Reynolds numberin reality, and in our scaled analogue experiments outlined below, the apparent viscosity of the magma varies spatially, and the flow profile is not parabolic. Next, we must choose a HEC solution to scale to the magma. Let us try X c = 0.8 wt% at 20°C. We use the Eqs. (13), (14), and (15), along with the appropriate parameter values, given in Section 3.5.1, to determine the rheological properties of the solution: η 0 , c, and p. We determine μ n = η 0 /(1 − ϕ) −1 hence η n0 (from Eq. (24)) and fit the dimensionless Cross model for HEC (Eq. (26)) to the dimensionless Cross model for the bubbly magma (Eq. (22)). In this case we find a good fit for scale factor b = 2.6 (maximum mismatch between the dimensionless rheological models is b15% for Ca ≤ 10) with viscosity ratio μ r = 290.0, and strain-rate factor bc/K = 2.2.
We now check that the Reynolds number is appropriately scaled, following the same approach as we did for the bubbly magma, but using the values appropriate for the scaled analogue fluid and experimental kithere we assume that the experimental analogue dyke has width 5 cm. The strain-rate at the wall of the analogue dyke is found from the strain-rate factor: γ • ¼ Ca=1:5 ¼ 4:6 (since Ca = 10 at the wall). From this, we can calculate the pressure gradient, flux, velocity, hence finally the Reynolds number. In this case, Re = 1.3. Comparing this with the value computed for the bubbly magma (Re = 1.3), we conclude that the experiment is well scaled for dynamic regime. This solution is well scaled for both the shear thinning rheology, and the dynamic Fig. 9. (a) Dimensionless comparison of the rheologies of a bubble suspension with bubble fraction ϕ = 0.4 and polydispersivity index m = 1.2, and a HEC solution with concentration X c = 1 wt% at temperature T = 20°C. The viscosity of the HEC is normalized so that η n0 = η r0 (Eq. (26)). (b) Cross number is scaled to improve agreement between the two curves for Ca b 10 (Eq. (31)). regime; hence, it is suitable for use in the proposed analogue experiments, and we expect the spatial gradient in viscosity and the shape of the velocity profile across the flow to match the natural case. Note that it is not always possible to scale both rheology and dynamic regime for a given magma and HEC concentration. For instance, in the current example a solution of X c = 1.5 wt% does not scale well for reasonable experimental temperatures.

Potential range of application
HEC is shear-thinning over a wide range of strain-rates. This means that it has potential for use as a shear-thinning analogue in experiments addressing a range of different flow problems in volcanology. To illustrate this, we compare our results to previous analogue studies in the literature that, together, span a wide range of strain-rates; these include: (1) lava flow characteristics (Soule and Cashman, 2005); (2) the ascent of gas slugs in a conduit (Llewellin et al., 2011); and (3) decompression driven fragmentation and permeable outgassing within a shock tube (Namiki and Manga, 2008;Namiki and Manga, 2005). For each of these studies we calculated the typical range of strain-rates encountered and plotted the ranges on the viscosity-strain-rate curves for all HEC concentrations at 20°C (Fig. 10). These ranges are all at higher shear-rates than the zero-shear viscosity plateau, meaning that the fluid viscosity will decrease with increasing shear. Therefore, if existing experimental methods were used with HEC the effect of a non-Newtonian shear-thinning rheology could be investigated.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  (Soule and Cashman, 2005); the ascent of gas slugs in a conduit (Llewellin et al., 2011) and decompression driven vesiculation Manga, 2008, 2005) Note that these black lines denote shear-rate ranges only, not viscosity.