Journal of Geophysical Research: Solid Earth The Gravitational Stability of Lenses in Magma Mushes: Conﬁned Rayleigh-Taylor Instabilities

In the current paradigm, magma primarily exists in the crust as a crystalline mush containing distributed melt lenses. If a melt-rich (or ﬂuid) lens is less dense than the overlying mush, then Rayleigh-Taylor (RT) instabilities will develop and could evolve into spheroids of ascending melt. Due to contrasting melt-mush rheologies, the theoretical RT instability wavelength can be orders of magnitude larger than the magmatic system. We explored how this conﬁnement aﬀects the gravitational stability of melt lenses through laboratory experiments with pairs of liquids with one layer much thinner and up to 2 . 2 ⋅ 10 5 times less viscous than the other; we extended the viscosity ratio to 10 6 with linear stability analysis. We found the growth rate of a bounded RT instability is approximately Δ 𝜌 gD 6 𝜋𝜇 2 , where Δ 𝜌 is the diﬀerence in density between the ﬂuids, g is gravity, D is the container diameter, and 𝜇 2 is the viscosity of the thicker viscous layer. This diﬀers from the unbounded case, where the growth rate also depends on the thickness and viscosity of the thin, low-viscosity layer. Applying the results to melt lenses in magmatic mushes, we ﬁnd that for the ranges of expected rheologies, the timescales for development of the instability, and the volumes of packets of rising melt generated span very wide ranges. They are comparable with the frequencies and sizes of volcanic eruptions and episodes of unrest and so suggest that RT instabilities in mush systems can cause episodic volcanism.


Introduction
A major challenge of modern volcanology concerns subsurface magma transport and accumulation. Conceptual models are emerging that depict subsurface systems as large uneruptible crystalline networks (mushes) containing heterogeneously distributed pockets of eruptible magma and exsolved volatiles that can extend deep in the crust and down to the mantle (e.g., Bachmann & Huber, 2016;Cashman et al., 2017). The dynamics of igneous mush systems has become a dominant theme in contemporary magma physics (e.g., Bergantz et al., 2017;Parmigiani et al., 2014) and a key feature of interpretations of geophysical, geochemical, and petrological data (e.g., Jaxybulatov et al., 2014;Putirka, 2017).
One aspect of igneous mush dynamics is the development of buoyancy instabilities related to intrusion of new magma or segregation of melt-or fluid-rich layers within a mush, leading to magma or fluid ascent through the mush. Magmas, melts, and fluids are commonly less dense than the overlying mush; therefore, Rayleigh-Taylor (RT) instabilities develop naturally wherever buoyant layers form. However, for some conditions (e.g., sufficiently high mush viscosity) the growth rate may be sufficiently slow that other processes (e.g., solidification due to cooling) dominate.
The viscosity contrast between a buoyant layer and an igneous mush is typically very large. For example, if we consider mushes with effective viscosities from 10 13 Pa⋅s for melt-rich mush (∼40% melt) to 10 17 Pa⋅s for melt-poor mush (<10% melt, Costa et al., 2009;Lejeune & Richet, 1995), and melt (magma) lenses with viscosities of 1 − 10 5 Pa⋅s (basalt to wet rhyolite), then the viscosity ratios are of order 10 8 -10 17 . There is a theory for RT instabilities for infinite horizontal layers with very high viscosity ratios (Whitehead & Luther, 1975), which has been verified in experiments with viscosity ratios of up to a few hundred. Although never tested experimentally, this theory should be applicable to the much higher viscosity ratios expected in lens-mush magmatic systems. However, a consequence of the high viscosity ratios is that the horizontal dimension of a magma reservoir is commonly much less than the theoretical fastest-growing wavelength assuming an infinite horizontal layer ( Figure 1). Thus, it is necessary to consider cases in which the buoyant layer is confined by boundaries separated by lengths much less than the optimum wavelength. The purpose of this study Figure 1. Sketch of an eruptible melt-rich layer within a much more viscous crystal mush. Rayleigh-Taylor instabilities arise due to density difference; however, the theoretical wavelength of instability c may be much larger than the melt layer diameter D.
is to address this issue through experimental and theoretical approaches and then to apply the results to understanding RT instabilities in igneous mushes. The suggestion that RT instabilities can control the frequency of volcanic eruptions is explored.

RT Instabilities
A RT instability occurs when a dense fluid overlies a less dense fluid. Such instabilities can arise in magmatic systems when a layer of buoyant melt is trapped within a denser crystalline mush. The full dynamics of RT instability are too complex to be described analytically. Nonetheless, linear stability analysis (LSA) provides good insight into the early stage, during which small initial perturbations of the interface grow exponentially (Waddell et al., 2001).
Mush-melt RT instabilities involve very high viscosity ratios. The relevant low Reynolds number RT instability formalism for such high viscosity contrasts was developed by Whitehead and Luther (1975). They analyzed the case of a thin layer of light Newtonian fluid (density 1 , viscosity 1 ) below a deep layer of denser Newtonian fluid ( 2 , 2 ); both fluids are horizontally infinite ( Figure 1). In this scenario, for a free slip boundary condition (BC) below the thin layer, and a large viscosity ratio ( = 2 1 ≫ 1), there exists a fastest-growing mode, the critical mode, with wavelength c and growth rate n c given by where h is the thickness of the thin layer and Δ = | 2 − 1 | is the density difference. It is also assumed that the wave-like perturbations are small, with amplitude less than ∼ 0.4 c . These equations have been modified to investigate the thickening of a buoyant layer at a fixed rate (i.e., constant dh dt ; de Bremond d' Ars et al., 1995). The theory has been tested and verified experimentally and numerically for viscosity ratios up to = (10 2 ) (de Bremond d' Ars et al., 1995;Whitehead & Luther, 1975).
Equations (1) and (2) assume that both layers are horizontally infinite, but for a finite domain, for sufficiently large h and , the theoretical fastest-growing wavelength, c , will be greater than the extent of the layers and so is not physically possible. Whether the domain is finite or infinite, all perturbations of all wavelengths will grow but the instability will develop with the wavelength that grows fastest. So we need to consider the growth rate of all modes with wavelengths equal to or smaller than the horizontal extent of the layers. In general the growth rate of a mode with wavelength and dimensionless wave number K = 4 h is (Whitehead & Luther, 1975): .
(3) as a function of dimensionless wave number K, for = 10 6 . Any mode with K < K D (gray area) cannot form because its wavelength is greater than the diameter of the container. As an example we illustrate this for K D = 0.13, which corresponds to h D = 0.01. For this case, the critical mode (K c ) cannot form. For large K,

Journal of Geophysical Research: Solid Earth
The normalized function n(K) is plotted in Figure 2 for = 10 6 . It reaches a maximum at K c = 4 h c , such that n(K c ) = n c (equation (2)). For K ≫ K c and ≫ 1, the growth rate decreases with increasing K, according to If the RT instability occurs in a container with horizontal dimension D < c , the tank walls will prevent the growth of all modes with K < K D = 4 h D . Graphically, any mode within the gray shaded area in Figure 2 is inaccessible. Thus, in the presence of lateral boundaries with K D > K c (i.e., c > D), the fastest allowed growth rate is n(K D ). It is useful to define a confinement parameter = c D , such that the system is confined for > 1. Using the definitions of and c , we can rewrite K D = 2.88 −1∕3 , therefore yielding This suggests that, in the case of a laterally confined RT instability, the growth rate is independent of both the thin layer thickness h and the viscosity ratio , though equation (5) requires experimental validation.

Are Mush-Melt Systems Confined?
The application of equations (1) and (2) to magmatic systems requires that melt layers can be considered horizontally infinite. To assess the condition for which this assumption is valid, we compute the confinement parameter = c D for natural magmatic systems. A system will be confined (or "clipped"; Burgisser & Bergantz, 2011) for > 1, but it can be considered unconfined when ≤ 1.
We examine the confinement of three example scenarios: (a) a silicic, (b) an andesitic, and (c) a basaltic system, with viscosities and densities as indicated in Table 1. For each of the three systems, we calculate for a range of melt layer widths (diameters if a circular lens) 500 ≤ D ≤ 2⋅10 4 m and heights 0.1 ≤ h ≤ 1, 000 m (Figure 1), which encompass typical values assumed in theoretical models (e.g., Annen et al., 2015;Bachmann & Bergantz, 2004) and consistent with estimates from geophysical surveys of natural systems (e.g., Lees, 1992;Tarasewicz et al., 2012). The majority of mush-melt systems are laterally confined (Figure 3). In fact, for the rheologies explored (Table 1), even an extremely wide (D = 20 km) layer of silicic (viscous) melt would be confined if thicker than 4 m. This result suggests that we should not simply apply equations (1) and (2) to typical lenses in magma mushes because the mechanical constraints from the side boundaries will prevent the growth of the critical mode with wavelength c . Instead, we expect the dynamics of the confined RT instability to be given by the fastest mode whose wavelength fits within the layer width. Therefore, we need a modified theory for RT instability that accounts for wall effects.

Fluids
To model mush-melt RT instabilities, we perform high viscosity contrast and laterally bounded experiments in cylindrical tanks, at low Reynolds number. We use glucose syrup (42DE-GL0106 from Ragus) as our viscous fluid and aqueous potassium carbonate (K 2 CO 3 ) solutions as our low-viscosity fluids. Physical properties of these two fluids are summarized in Table 2 at both 18 and 8 ∘ C. Both fluids are Newtonian, and viscosity ratios in our experiments are in the interval 4.7 ⋅ 10 4 ≤ ≤ 2.2 ⋅ 10 5 , a range which extends previous experimental data to much higher ratios.
Viscosity measurements were performed on a HAAKE RheoStress 1 rheometer (Thermo Fisher Scientific), with a concentric cylinders sensor system. Glucose viscosity at 8 ∘ C was too high to be measured directly as the rheometer was not calibrated to such high values. Instead, we measured viscosity from 18 to 28 ∘ C, and, assuming an Arrhenius model, we fitted an exponential curve through these data and extrapolated it down to 8 ∘ C. We measured glucose density via Archimedes' principle, by weighing a calibrated 10-cm 3 sinker (Mettler Toledo) in both air and glucose. The very low viscosity of the K 2 CO 3 solutions allowed us to measure density by weighing 200 ml of solution in a volumetric flask.

Experimental Methods
We use two tank sizes to investigate the effect of lateral confinement. All tanks are made of transparent Perspex and are cylindrical to avoid corner effects.

Narrow Tank Experiments
The narrowest tanks have a diameter D = 5.00 ± 0.02 cm and height H = 40.00 ± 0.05 cm ( Figure 4). Three tanks of this geometry were made in order to efficiently run concurrent experiments. The confinement factors ( = c D ) for experiments in these tanks are 10 ≤ ≤ 210; hence, wall effects are likely to be very important in all experiments with D = 5 cm.
We conducted experiments in two temperature-controlled rooms, at 18 ∘ C (laboratory) and 8 ∘ C (cold room), to sample a wider range of viscosity ratios. We first pour more than 30 cm of glucose syrup into the tank, cover the top of the tank, and leave it to rest until all the air bubbles have escaped. The uppermost part of the syrup usually dehydrates and stiffens. We thus remove this stiff layer a few minutes before starting the experiment. This is long enough for the disturbed interface to flatten but short enough to avoid significant drying. The experiment is started by delicately pouring dyed K 2 CO 3 solution on top of the glucose layer. Pouring takes a few seconds, which is much faster than the time for an RT instability to develop (always >3 min). This configuration (viscous fluid underneath) is flipped compared to magmatic systems with a buoyant lens under a more viscous mush. Nonetheless, the dynamics of the instability will be identical because the driving force, buoyancy, is independent of which fluid is on top of the other; only the density difference Δ is important (Whitehead & Luther, 1975). All experiments are recorded with a fixed camera.
In most experiments, there is air directly above the K 2 CO 3 solution, forming a free-slip BC above. To check the influence of the BC, we repeated some experiments with a no-slip upper BC by placing a circular perspex lid on top of the K 2 CO 3 layer just after pouring it. The lid diameter is 0.4 mm less than the tank to allow air escape during placement.

Wide Tank Experiments
We perform similar experiments in a wide tank with diameter D = 28.70 ± 0.05 cm and height H = 39.20 ± 0.05 cm, achieving 1.2 ≤ ≤ 6.9. Although > 1, we expect wall effects to be small and our results to approach the theoretical predictions of Whitehead and Luther (1975). The experimental procedure is identical to the narrow tank case, with the exception that we pour the K 2 CO 3 solution through a nozzle with seventy-two 1.5-mm-diameter holes to reduce disturbance of the interface. Again, we conduct experiments at both 18 and 8 ∘ C. We, however, did not run any no-slip upper BC experiment with the wide tank.

Linear Stability Analysis
To complement our experimental investigation, we perform an linear stability analysis (LSA), based on the work of Sweeney et al. (2013). They theoretically investigated RT instabilities in a narrow, finite cylinder (no-slip BC), but their published results cannot be directly applied to our scenario for two reasons. First, they only examined cases where the two fluid layers have equal thicknesses, and second the viscosity ratios considered are too low ( ≤ 10 2 ). We use their numerical routines to reproduce our experimental conditions and to extend results to higher values than achieved in the laboratory. Full technical details can be found in Sweeney et al. (2013). Note. All symbols are defined in the text. Unless indicated by "Wide Tank", all experiments were performed in a narrow tank with D = 5.00 ± 0.02 cm. The final column indicates n c n obs for n c calculated by equation (2) and n obs determined by linear stability analysis.

Results
We conducted a total of 40 experiments, sampling the range 1.2 ≤ ≤ 210. The viscosity ratios (4.7 ⋅ 10 4 ≤ ≤ 2.2 ⋅ 10 5 ) are more than 2 orders of magnitude higher than previous large viscosity contrast RT experiments (e.g., de Bremond d'Ars et al., 1995;Huppert et al., 1984). A summary of the experimental parameters and corresponding results is presented in Table 3. We performed an LSA for each set of experimental parameters (Table 3). To complete and extend our data set, we also carried out LSA in the ranges 2.3 ≤ ≤ 7.0 and 500 ≤ ≤ 10, 000, with sets of parameters that have not been or could not be experimentally investigated. The conditions and results for these additional LSA runs are summarized in Table 4. We thus obtain theoretical growth rates for the confined RT instabilities considered, which can be compared to our experimental observations and to unconfined theory. Figure 5 shows the time series of two typical experiments in a narrow tank: one at room temperature with air above (Figure 5a) and one at 8 ∘ C with a no-slip upper BC (Figure 5b). Regardless of the tank size or the BC, the overall dynamics are similar. First, a single, small protrusion of K 2 CO 3 solution forms at the interface between the two fluids. The protrusion grows into an ellipsoidal to spheroidal pocket of solution. When all of the initial Figure 5. Time series of the initial stage of (a) experiment 25, at 18 ∘ C with h = 1 cm, Δ = 50 kg/m 3 , = 6.6 ⋅ 10 4 , and = 35 and (b) experiment 40, at 8 ∘ C with h = 1.7 cm, Δ = 52 kg/m 3 , = 2.2 ⋅ 10 5 , = 89 and a no-slip boundary condition. The initial interfaces are highlighted and an example of how amplitude is measured is provided. The location of the initial protrusion is sensitive to asperities on the syrup surface or disturbance induced by pouring the K 2 CO 3 solution; hence, the protrusion is sometimes initially slightly off centered (Figure 5b for instance). Replicate experiments show that these deviations do not significantly affect the results: The growth rate of the protrusion is similar, and the K 2 CO 3 solution layer always becomes thickest near the center of the tank well before it evolves into a spheroid sinking through the syrup.
The amplitude of the perturbation, defined as the vertical distance from the initial liquid-liquid interface location to the bottom of the protrusion (see Figure 5a), was measured as a function of time for each experiment ( Figure 6). For all experiments there is an initial exponential increase in amplitude before transitioning to linear growth. This transition corresponds to the stage at which linearized stability theory (e.g., Whitehead & Luther, 1975) is no longer valid. We therefore obtain an exponential fit of the form y(t) = Ae n obs t for the initial growth of each experiment. The n obs parameter in the exponential then corresponds to the growth rate of our experimental RT instability and can be compared to theoretical unconfined predictions.  In summary, our results show that bounded growth rates (n obs ) are reduced compared to the unbounded theoretical values (n c ) by approximately a factor of . For simplicity, in our applications to magmatic system we take n c n obs = , which is within 10% of the lines of best fit to our data ( Figure 7) and introduces negligible error compared to the uncertainties on the magma and mush properties that affect RT instabilities.

Confined RT Instability Dynamics
Our experimental results with 1.2 ≤ ≤ 210 and 4.7 ⋅ 10 4 ≤ ≤ 2.2 ⋅ 10 5 demonstrate that the dynamics of a confined RT instability is qualitatively similar to the unconfined case. Indeed, the instability starts with an exponential growth, followed by a linear growth (e.g., Figure 6), a phenomenon well documented for the unbounded case (e.g., Waddell et al., 2001). Moreover, the shape of the perturbations is similar to previous experimental studies (e.g., Waddell et al., 2001;Whitehead & Luther, 1975;Wilkinson & Jacobs, 2007). Our instabilities differ from the unbounded case because only a single protrusion forms, as opposed to multiple, uniformly spaced bulges, as observed by Whitehead and Luther (1975) and de Bremond d'Ars et al. (1995) for instance.
Our experimental and analytical results indicate that as approaches unity, there is a transition in the instability growth rate from n c n obs = when the system is bounded towards n c = n obs , characteristic of an unbounded ( ≪ 1) system (i.e., n obs matches the theoretically predicted value for a laterally infinite system). Based on these results, we estimate a characteristic RT instability timescale, that is, the time for the instability amplitude to increase by a factor of e (≈ 2.72), by In contrast to the unbounded case, the confined timescale does not depend on either the initial layer thickness h or the viscosity ratio because the instability wavelength is set by D. Rather, the controlling parameters are the diameter of the layer D, the viscosity of the upper layer 2 , and the density contrast between the fluids Δ . Comparing equations (5) and (6), we notice that our experimental growth rates are a factor of 3/2 slower than estimated with the linear stability theory of Whitehead and Luther (1975) for = D. This is most likely because the tank walls add a no-slip BC (i.e., will exert significant drag on the fluids) that is not accounted for in equation (5).

Application to Igneous Systems
There are multiple mechanisms where buoyant layers can form within igneous mush systems. One mechanism is by replenishment with new magma spreading out within or at the base of a mush. Initially, the new magma layer may intrude as a denser layer at the base of the reservoir or at its neutral buoyancy level within the reservoir. In either case the magma layer may become less dense with time by, for example, crystallization and differentiation with dense components segregating to the base of the flow. Volatile exsolution can also increase the buoyancy of the layer (Huppert et al., 1982). Additionally, the overlying mush may be heated from below (e.g., Burgisser & Bergantz, 2011;Couch et al., 2001) or be fluxed by volatiles released from the replenishing magma layer (Bachmann & Bergantz, 2006) to develop a zone of reduced density. Another mechanism of layer formation is through dynamic melt percolation due to mush compaction, which can result in formation of melt-rich regions. One-dimensional models of porous media flow and compaction predict the development of multiple melt-rich regions (e.g., Jackson et al., 2003;Solano et al., 2012). Intrinsically, 1-D models cannot include RT instabilities of growing melt layers. Likewise, exsolved magmatic volatiles can migrate through mushes and accumulate as fluid layers (Christopher et al., 2015). Each of the above scenarios for melt-or fluid-layer generation could be modeled individually and in detail, but this is beyond the scope of the current paper. Here we investigate the stability of buoyant magma layers of fixed thickness or growing at a fixed rate to provide a first-order understanding of the likely timescales and length scales that can be expected. Our analysis is restricted to RT instablities; we have not explored the potential importance of melt-mush flow mechanisms and smaller-scale instabilities related to the porous or brittle nature of the mush, which could develop along with or instead of RT instabilities (e.g., Connolly & Podladchikov, 2013;Oppenheimer et al., 2015;Sandnes et al., 2011;Schmeling et al., 2017;Scott & Stevenson, 1984, 1986. Figure 8 presents estimates of confined RT instability timescales for a range of mush viscosities 2 and melt layer diameters (or widths if not circular) D relevant to igneous systems. We use the instability timescale defined in equation (6), with calculated as c ∕D using c defined in equation (1), and physical parameter values listed in Table 1. The wide range of calculated timescales (10 1 -10 4 years) match the variability of timescales observed for volcanic processes, from small and frequent events to large and rare events (e.g., caldera forming eruptions). On the other hand, the very long timescales involved with the slowest instabilities imply that other processes (e.g., solidifaction by cooling) could occur faster and inhibit the onset of RT instability.

Comparison to Cooling Timescale
First, we consider a case where there is a strong temperature difference between the melt and the mush and compare our RT instability timescale to a characteristic cooling timescale. Such a scenario could occur when a cold mush system is replenished by hot magma from depth. In this case, if the cooling timescale is shorter than the RT instability timescale, the melt layer will freeze before an instability can develop. For simplicity, we assume cooling occurs via conduction only and define a characteristic cooling timescale cool = h 2 , where h is the melt layer thickness and is the thermal diffusivity. We use = 6 ⋅ 10 −7 m 2 /s (Annen et al., 2006;Romine et al., 2012;Whittington et al., 2009). Figure 9 shows the lines where cool = RTI as h and D are varied and for two mush viscosities. Above these lines, RT instabilities have time to develop, whereas under the lines, the melt layer will freeze before the instability develops significantly. The results are dramatically different depending on which mush viscosity is considered. For a melt-rich mush ( 2 = 10 13 Pa⋅s), RT instabilities develop sufficiently quickly that the required thickness for the instability to develop faster than conductive cooling is only 1.4-8.9 m, depending upon the layer diameter. This range, however, becomes 140-900 m for a near-solidus mush ( 2 = 10 17 Pa⋅s). The mush rheology therefore exerts a critical control on the necessary thickness to develop instability and hence the volume of eruptible material ascending through the mush as a result of the instability.
This simple analysis omits two potentially important but opposing mechanisms: (a) we only consider cooling via conduction and ignore convection, which can speed up cooling, and (b) the hot intrusion can reheat the mush, thus decreasing its viscosity and speeding up RT development. To get some insight about when convection may occur, we calculate the Rayleigh number Ra = g 1 ΔTh 3 1 , where is the thermal expansion of the melt and ΔT is the temperature difference between the melt and the mush. We use = 5 ⋅ 10 −5 K −1 and ΔT = 100 K and densities and viscosities for the three compositions considered in Table 1 (  . Isolines cool = RTI as a function of melt layer width and thickness and for two mush viscosities. Above these lines, there is sufficient time for the Rayleigh-Taylor instability (RTI) to grow, whereas under the lines, the melt layer will freeze before the instability develops significantly.
as the critical value for the onset of convection (Turcotte & Schubert, 1982), this calculation suggests that all the layers previously considered would convect, except for extremely thin silicic layers (h < 49 cm ). Therefore, the cooling timescale could be faster than assumed in Figure 9, which would raise the cool = RTI lines to greater h for a given D. Concerning reheating of the mush, Burgisser and Bergantz (2011) showed that RT instabilities could lead to significant mush overturn ("unzipping"), yet their analysis is restricted to low mush viscosities ( 2 = 10 6 − 10 12 Pa⋅s) and considers high temperature differences (ΔT = 45-500 K). Therefore, the importance of this mechanism has yet to be tested for our scenario.

Comparison to Accumulation Timescale
The previous discussion assumes that melt accumulation into a layer occurs much faster than the instability timescale, such that a melt lens of any thickness is a feasible starting condition for considering RT instabilities. Next, we relax this assumption and explore the relation between melt accumulation and RT instability timescales.
De Bremond d' Ars et al. (1995) investigated the case of a horizontally extensive buoyant layer that thickens and a constant rate,ḣ = dh dt . They showed that blobs of buoyant fluid form and rise away from the base when the layer growth ratėh h and the instability growth rate 1 RTI are equal. For the formation of a new layer (i.e., starting with h = 0), the layer growth ratėh h is initially very large and then monotonically decreases as h increases. Thus, for any given melt input rateḣ, there exists a maximum layer thickness beyond which melt will be removed faster via RT instabilities than it is added. We can estimate this thickness, the corresponding magma volume, and the time required to assemble it for a range of parameters. First, we calculate the time at whicḣh h = 1 RTI . We can then compute the layer thickness using the constant melt input rate and corresponding volume, assuming a cylindrical shape. We use the same criterion as de Bremond d' Ars et al. (1995) but account for lateral confinement when relevant. Figure 10 shows the accumulation time and associated melt volume for a silicic melt under a melt-rich mush ( 2 = 10 13 Pa⋅s; Figure 10a) and a near-solidus mush ( 2 = 10 17 Pa⋅s; Figure 10b), as a function of the layer diameter D and the melt input rateḣ. The chosen range of input ratesḣ span values from magma accumulation models and field studies (e.g., Karakas et al., 2017;White et al., 2006). The kinks in the lines for melt input rates ofḣ = 0.1, 1, 10 m/year in Figure 10a correspond to the transition from confined to unconfined instability regime with increasing D. For conditions where the RT instability is not confined (i.e., < D), the accumulation times are independent of D but depend onḣ (forming horizontal lines in the upper panel of Figure 10a). In confined scenarios (some combinations ofḣ and D in Figure 10a and all conditions plotted in Figure 10b), the accumulation times are independent ofḣ, hence the collapse into a single line for a given 2 . The choice of mush viscosity has a major impact on the results and allows us to recover a large range of timescales (0.3-4 ⋅ 10 4 years). The corresponding volumes range from 10 5 up to 10 12 m 3 .
Our calculations suggest that RT instabilities could play an important role in controlling the size and frequency of volcanic events. A working hypothesis is that volcanic eruptions and episodes of volcanic unrest are the consequences of these instabilities. The calculated timescales and volumes are comparable with natural volcanic values (e.g., Pyle, 2015) and span the full range between small mafic eruptions (e.g., Strombolian type) to the largest magnitude eruptions (e.g., caldera-forming eruptions and flood basalts). The time needed for an instability to develop allows melt to accumulate in large layers and corresponds to a dormant period (e.g., Rougier et al., 2018;Sheldrake et al., 2016), whereas the instability may destabilize the system and produce a period of unrest or an eruption. Additionally, successive instabilities without eruption could yield larger accumulated volumes (e.g., .

Crystal Mush Rheology
A limitation in the application of our experimental and theoretical results to natural systems is the assumption of Newtonian mush rheology. Our experiments involved Newtonian fluids, but natural crystal mushes are generally shear-thinning and may have a yield strength, that is, a minimum stress required for flow to occur (e.g., Hoover et al., 2001;Kerr & Lister, 1991;Saar et al., 2001). With Newtonian fluids, an RT instability will always develop in response to a denser fluid overlying a less dense fluid, although in some scenarios it will grow slowly enough that it can be neglected over the timescale of interest. The same is true of shear-thinning fluids; here the relevant viscosity for the early stage of the instability is the (high) viscosity in the limit of the shear rate approaching 0. With a yield strength, however, the buoyant force from the melt layer Δ gh has to exceed the mush yield strength 0 ; otherwise, the mush will act as a solid (e.g., elastic) body, preventing the growth of the RT instability, akin to initiation of thermal convection of a fluid with a yield strength (Balmforth & Rust, 2009). True yield strengths (i.e., a minimum stress for any flow to occur) may not exist (e.g., Barnes, 1999); however, if the effective viscosity of a mush at low stresses is extremely large such that it has an apparent yield strength, then the instability will grow so slowly as to be negligible unless Δ gh > 0 .
We can estimate the minimum thickness required for an RT instability to develop when the mush has a yield strength as h min = 0 Δ g . For mush yield strengths in the range 0 = 10 5 -10 6 Pa (Castruccio et al., 2013;Lejeune & Richet, 1995) and Δ = 300 kg/m 3 , we find h min = 34-340 m. This crude estimate suggests that the required thickness to overcome a yield strength could be of the same order of magnitude as typical melt layer thicknesses. Mush strength could thus facilitate the accumulation of melt lenses by impeding RT instability development for the thinnest layers. This effect could be enhanced if crystals have grown together and bonded in a stagnant mush, producing a yield strength greater than 10 6 Pa.

Conclusion
Rayleigh-Taylor instabilities occur naturally in magmatic systems when buoyant melt (or magmatic volatile phase) is trapped underneath a denser crystalline mush. For a wide range of expected viscosities, the large viscosity contrast between the mush and melt lens means that the theoretical fastest-growing wavelength is unfeasibly large and so the wavelength of the instability is the largest available: the diameter of the lens. This lateral confinement means that the growth rate of the instability is reduced compared to the theoretical unconfined scenario. Importantly, if confined, the instability growth rate no longer depends on the thickness or viscosity of the lens; rather, it depends only on the diameter of the lens, the viscosity of the mush and the density contrast between the lens and mush. The thickness of the lens will, however, play a role if the mush has an apparent yield strength because the instability will only initiate if the buoyancy stress, which is proportional to the lens thickness, is sufficient to overcome the yield strength. Thickness also matters in scenarios where the melt is hotter than the mush, as thinner lenses will cool more quickly, allowing less time for RT instability development before the magma is too crystalline to flow.
A fully developed confined RT instability transforms the melt lens into a spheroid of melt, which ascends through the mush. A lens of melt can only exist if it thickens faster than melt is removed by a RT instability. So both the timescale for transforming a lens into a rising spheroid and the volume of the spheroid will depend on the rate of input of melt into the lens as well as its diameter, the viscosity of the mush, and the density contrast between the lens and mush. We postulate that RT instabilities may play a role in regulating the size and frequency of volcanic eruptions and volcanic unrest. Using feasible ranges of the relevant parameters, we calculate timescales and volumes that span small, frequent mafic eruptions to the largest and much rarer caldera-forming and flood basalt eruptions. Also, multiple episodes of layer instability without eruption can lead to accumulations of larger magma volumes and provide one explanation of volcanic unrest. In a magmatic system composed of a vertically extensive mush containing multiple melt lenses, can the rise of a blob of melt or volatiles trigger a feedback to cause large-scale destabilitization? The growth rate of a laterally confined RT instability in a Newtonian mush does not depend on the thickness of the melt lens. So the ascent of a blob of melt into a more shallow lens will not cause an accelerated destabilization of that lens unless: (a) it causes the diameter of the lens to increase (and so increases the wavelength of the RT instability) or (b) it increases the thickness of the lens such that its buoyancy stress overcomes the yield strength of the mush above it. Mush rheology therefore is important for controlling where melt (and magmatic volatile phases) accumulate in layers, how much melt accumulates, and how frequently packets of melt are released.