Surface uplift is a phenomenon observed at man y acti ve volcanoes around the world (e.g. Biggs & Pritchard 2017 ). It may be due to magma movement underground, including episodes of renewed magma supply or recharge to subsurface reservoirs (e.g. Delaney & McTigue 1994 ;Dvorak & Dzurisin 1997 ;Dzurisin et al. 2006 ). Although we do not understand these processes with certainty, we can use geodetic observations to make inferences about the state of magma storage reservoirs as well as associated magma supply or recharge. The nature of this deeper magma supply process, ho wever , often remains obscure, non-unique and therefore treated in a highly idealized manner.
Our study focuses on inflationary surface deformation during the past decade at Cord ón Caulle volcano in Chile (Delgado et al. 2016(Delgado et al. , 2018 ; ;Delgado 2021 ). Cord ón Caulle is of particular interest because it has erupted crystal-poor rhyolite and contains evidence of an active crystal mush through the presence of crystal-rich mafic enclaves within the erupted rhyolite. The enclaves have a glassy matrix of identical composition to the er upted cr ystal-poor rhyolite (W inslo w et al. 2022 ). It has therefore been inferred that the erupted rhyolite represents melt that se gre gated from an underlying crystal mush, as is envisaged for the rhyolite-melt crystal mush model (Bachmann & Bergantz 2004 ;Hildreth 2004 ;Hildreth & Wilson 2007 ). Cord ón Caulle, whose eruptions are undoubtedly fed from a crystal-poor rhy olite reserv oir, represents a v olcanic system where the potential existence of a crystal mush can be tested against obser vations of inflationar y defor mation. Specifically, one can test whether the observed inflation could be the consequence of recharge of the crystal-poor rhyolite by melt se gre gation from an underlying crystal mush.
The rhyolite-melt crystal mush model is part of a conceptual model for lithospheric magmatism wherein asthenosphere derived Melt se gre gation from a crystal mush 611 basalt feeds transcrustal magmatic systems in which fractionation and hybridization processes e ventuall y produce rhyolite melts at shallow crustal levels (Hildreth 1981 ), in essence similar to present views of transcrustal magma transport and storage systems (e.g. Annen et al. 2006 ;Bachmann & Huber 2016 ;Cashman et al. 2017 ;Sparks et al. 2019 ). Although crystal-melt fractionation is believed to take place over a wide range of depths, it is thought that the temporal persistence and growth of such systems leads to ele v ated upper-crustal geothermal gradients that favour the formation of shallow bodies of high-silica rhyolite-melt by step-wise extraction by porous flow from underlying crystal mushes (Hildreth 2004 ;Hildreth & Wilson 2007 ;Huber et al. 2019 ).
No evidence for the existence of such a transcrustal system currently exists at Cord ón Caulle and considering how primitiv e the enclav es from the 2011-2012 eruption are, fractionation processes at Cord ón Caulle may be restricted to relati vel y shallo w depths (W inslo w et al. 2022 ). Nev ertheless, the enclav es do point to the existence of a crystal mush, perhaps at relati vel y shallow depths and limited v ertical e xtent. Here we assess whether the inflation at Cord ón Caulle during the past decade can be reconciled with melt se gre gation from such a crystal mush body.
We do so using a numerical model for magma transport and storage that is based on the framework of poroelasticity (e.g. Biot 1941 ;Rice & Cleary 1976 ;Wang 2000 ;Cheng 2016 ), which has been proposed to describe the mechanics occurring within a magma mush over shorter timescales than viscous or plastic deformation processes associated with compaction (Gudmundsson 2012 ;Liao et al. 2018 ;Mullet & Segall 2022 ). Hereby the crystalrich mush is treated as a porous, elastically deforming medium with interstitial melt. The model simulates how the decrease in pressure of the crystal-poor rhyolite reservoir, caused by eruptive magma withdrawal, led to upward flow of interstitial melt within the underlying crystal mush, resulting in recharge of the crystalpoor rhyolite and thus the observed inflation. We use the model to estimate the parameters under which it can reproduce the observed surface deformation and assess the extent to which they support the existence of a magma mush system beneath Cord ón Caulle.
Cord ón Caulle is part of the Puyehue-Cord ón Caulle volcanic complex (PCC), located within the Andean Southern Volcanic Zone (Fig. 1 (a); Singer et al. 2008 ). PCC is a long lived system which has produced ≈7 km 3 of erupted material in the past 16.5 ka. Cord ón Caulle is situated within a NW trending graben in the complex, and has produced three eruptions during the past century, in 1921-1922 and 1960 with eruptiv e v ents along graben-bounding faults, and in 2011-2012 (Fig. 1 b). All three eruptions were rhyolitic to dacitic and homogeneous in composition (Castro et al. 2013 ;Singer et al. 2008 ), suggesting that they were tapping different parts of a compositionally zoned magma storage system that is either contiguous or compartmentalized with some degree of interconnection, as indicated by spatiotemporal patterns of surface deformation (Jay et al. 2014 ;Delgado 2021 ).
On 4 June 2011, Cord ón Caulle began erupting with a 3-4-d-long Plinian to subplinian phase, with ash reaching 9-12 km, and a volcanic e xplosivity inde x (VEI) of 4-5 (Jay et al. 2014 ;Bonadonna et al. 2015b ). During the following week the ash plume decreased to 4-9 km in height (Bonadonna et al. 2015b ), and on June 15 crystalpoor rhyolitic lava began extruding from the vent, along with continued ash venting (Crozier et al. 2022 ). Seismic tremors ceased, possibly signalling the end of the eruption on 15 March 2012 (Tuffen et al. 2013 ). The exact end date is, ho wever , obscured since the lava flow field continued to evolve for months (Bertin et al. 2012(Bertin et al. , 2015 ; ;Tuffen et al. 2013 ;Coppola et al. 2017 ). Thus, the eruption may have ef fecti vel y ended b y Ma y 2012, when estimated la va discharge rates declined to < 5 m 3 s -1 (Coppola et al. 2017 ). In total, approximately 1.5 km 3 dense rock equivalent (DRE) of tephra and lava were e xtruded ov er the course of the eruption, with nearly 1 km 3 erupting e xplosiv ely within the first few days, and the lava flow making up the remaining ∼0.5 km 3 (Pistolesi et al. 2015 ;Castro et al. 2016 ;Delgado et al. 2019 ;Delgado 2021 ). Finally, a shallow ( < 200 m) laccolith of ∼0.8 km 3 also intruded syn-erupti vel y (Castro et al. 2016 ), bringing the total volume of magma withdrawn from the Cord ón Caulle magma storage reservoir to > 2 km 3 . This ultimately resulted in the deflationary deformation shown in Fig. 1 (c), which was preceded and followed by inflationary deformation, all of which will be discussed in more detail in Section 2.4.
The 2011-2012 lava is a crystal-poor rhyolite (Castro et al. 2013 ;Jay et al. 2014 ;Seropian et al. 2021 ) containing 1 volume per cent of mafic enclaves (W inslo w et al. 2022 ). These mafic enclaves have been interpreted to be pieces of crystal mush which contain interstitial rhyolite melt (W inslo w et al. 2022 ). The enclaves range in size from 2 to 20 cm and are comprised of 55-70 per cent interlocking crystals with the remainder as interstitial glass and vesicles (see Fig. 2 ). The enclaves' bulk composition is basaltic, due to an abundance of olivine, clinop yrox ene and plagioclase crystals, as well as minor amounts of Fe-Ti oxides. The interstitial glass of the enclaves has an identical composition to that of the crystal-poor host rhyolite. This fact, together with their crystallinity has led to the interpretation that most of the enclaves are pieces of a crystal mush that is spatially and genetically associated with the crystalpoor rhyolite body from which the 2011-2012 magma erupted. A possible scenario that would be consistent with the pre v alent view of how crystal-poor rhyolite is formed (Bachmann & Bergantz 2004 ;Hildreth 2004 ) would be a crystal mush overlain by crystalpoor rhyolite with the former somehow being incorporated into the rhyolite flowing towards the eruption-feeding dike or conduit (Hermes & Cornell 1981 ;W inslo w et al. 2022 ).
It should be noted that a small fraction of the enclaves have porphyritic textures, low crystallinity, plagioclase with disequilibrium textures, abundant microlites and a narrow range in mineral chemistry, all of which point towards these enclaves coming from direct mafic injection into the rhyolite layer (W inslo w et al. 2022 ). Their abundance is low, and it is unclear whether or to what extent their injection could be the cause of the inflationary deformation prior to the eruption and perhaps the eruption itself. Because the objective of our analysis is to test the possibility that the inflationar y defor mation following the er uption is due to melt seg regation from a crystal mush, we assume that this deformation is not caused (1921-1922), blue (1960) and green (2011-2012). Approximate locations of faults, based on Lara et al. ( 2006 ) and Delgado ( 2021 ) b y rene wed mafic injection. We ackno wledge, ho wever , that this possibility cannot be ruled out based on the available information.
Over the course of the eruption approximately 4.25 m of deflation was observed at Cord ón Caulle (Delgado 2021 ), using interferometric synthetic aperture radar (InSAR). Some spatiotemporal variability aside, deflation was essentially centred within the Cord ón Caulle graben (Jay et al. 2014 ;Delgado 2021 ;Novoa et al. 2022 ) and was consistent with magma withdrawal from a shallow magma body located within the graben (Delgado et al. 2019 ). Using a temporal sequence of interferograms, Delgado et al. ( 2019 ) found the line of sight (LOS) displacements near the centre of the graben. The resultant time-series of LOS displacements captures the temporal evolution of co-eruptive subsidence (blue diamond, Figs 1 b and c). Uncertainties exist due to lack of coherence between interferograms during the eruption and the spatiotemporal variability in the location of maximum deflation. They can, ho wever , still be used as a proxy for the relative rate of change in volume (pressure) within the magma reservoir.
Following the eruption there has been near continuous uplift at Cord ón Caulle, also approximately centred within the graben and with relati vel y limited spatiotemporal variability (Delgado et al. 2018 ;Delgado 2021 ). This uplift has been interpreted to be the consequence of magma recharge into the reservoir from below (Delgado et al. 2016(Delgado et al. , 2018 ; ;Novoa et al. 2022 ), although alternate interpretations may be feasible (e.g. Segall 2016 ). After analysing the post-erupti ve interfero grams, Delgado et al. ( 2016 ) produced a time-series of inflationary LOS displacements from a single location within the graben (red circle, Figs 1 b and c), which also can be used as a proxy for the relative rate of change in volume (pressure) within the magma reservoir. Although these data are given for a single location, the maximum deformation has shifted position around this average location. We also note that gradual inflation was occurring prior to the eruption at a slower rate and smaller magnitude than has been observed since the eruption (Delgado et al. 2018 ;Delgado 2021 ). The exact cause of this pre-eruptive inflation is unknown, but may have been related to deep-seated recharge. The LOS displacement data, recorded for a variety of satellite orbits, viewing geometries and slight spatiotemporal variations in maximum displacement, are scaled to provide a continuous timeseries as explained in Delgado ( 2021 ), to make the eruptive deflation continuous with the post-eruptive inflation (Fig. 1 ). For the purposes of our analysis, we treat the LOS displacements as a proxy for the rate of magma withdrawal and resultant change in pressure of the shallow magma reservoir that fed the 2011-2012 eruption. Although this proxy may not strictly describe the maximum deformation over this time span at Cord ón Caulle, it provides a sufficient proxy to compare against to test whether the inflationar y defor mation was caused by melt se gre gation from a crystal mush.
Our objective is to understand the inflationar y defor mation at Cord ón Caulle following the 2011-2012 eruption. In particular, we are interested in assessing whether the observed deformation could be due to melt se gre gation from a crystal mush, as it is envisioned within the the magma mush hypothesis for the formation of eruptible crystal-poor rhyolite (Bachmann & Bergantz 2004 ;Hildreth & Wilson 2007 ;Cashman et al. 2017 ). To this end, we developed a mathematical model that describes the poroelastic response of a crystal-rich magma mush to withdrawal of magma from an overlying layer of crystal-poor rhyolite that is contiguous with the mush. The system is presumed to exist within an elastically deformable host rock whose deformation is calculated using the conventional semi-analytical model for the opening of a sill-like dislocation (Okada 1992 ), which has been shown to provide a reasonab le appro ximation to deformation at Cord ón Caulle (Delgado et al. 2016 ). We emphasize that the conceptual model of our system comprises the magma mush plus an overlying sill of rhyolite melt. Rhyolite melt can be removed from this system due to eruption, which will result in a change in magma volume. We approximate the system as closed to magma recharge from below the mush. The aforementioned porphyritic enclaves, which are interpreted to be a consequence of injection of mafic magma, do allow for the possibility that inflation at Cord ón Caulle is due to the injection of magma from below the system. Although they are volumetrically insignificant in the erupted products, it is possible they make up a larger, trapped fraction in the mush. Moreover, the existence of the magma mush is of course contingent upon deep mafic recharge over time. The possibility of direct recharge from a deeper source through a dike or conduit has already been investigated by Delgado et al. ( 2016Delgado et al. ( , 2018 ) ), who found that an input of 0.15 km 3 into the sill through 2017 could produce the uplift. We therefore limit our analysis to the end-member case of recharge due to melt se gre gation from a magma mush in order to assess solely the extent to which the inflationar y defor mation during the past decade can be explained by this process.
We assume that the crystal mush acts as a porous medium, where the crystals constitute the porous matrix and the interstitial rhyolite melt is the pore fluid. The pressure within the overlying rhyolite acts as a pressure boundary condition on the mush. Pressure within the overlying rhyolite can change either due to magma withdrawal during eruption or due to magma addition from the underlying mush. A decrease in pressure of the overlying rhyolite, due to eruptive withdrawal, will induce pressure gradients in the interstitial melt that deviate from static, thereby inducing flow of interstitial melt within the mush. This results in a net flow of interstitial melt out of the mush into the overlying rhyolite.
If the resultant change in interstitial melt content were solely due to the compressibility of the interstitial melt, which would be the case for a completely rigid and undeformable crystal matrix, the volume of interstitial melt added or removed from the mush would be associated with no change in volume of the mush itself. On the other hand, if a change in interstitial melt content were solely due to the deformation of the crystal matrix, which would be the case if the interstitial melt is incompressible, then the volume of interstitial melt removed from the mush would be entirely compensated by a corresponding change in volume of the mush. In reality, any change in interstitial melt content is due to both melt compressibility and matrix deformation, where compressibility of the interstitial melt is due to exsolution and decompression of volatiles (Sparks & Cashman 2017 ). As a consequence, a redistribution of melt between mush and overlying rhyolite will result in a change in volume of the system. Moreover, the less deformable the crystal matrix, the larger the change in the volume of the system.
The theoretical framework of poroelasticity (e.g. Biot 1941 ;Rice & Cleary 1976 ;Wang 2000 ;Segall 2010 ;Cheng 2016 ;Liao et al. 2018 ) describes these principles. A change in volume of the mush therefore requires deformation of the crystal matrix, which is assumed to be instantaneous and fully recoverable. At the same time, the flow of interstitial melt follows Darcy's law and is timedependent. The flow can be described by a diffusion equation for pore pressure. As a consequence, an instantaneous change in pressure of the overlying rhyolite will result in the gradual diffusion of interstitial melt pressure throughout the mush, with any local volume change of the crystal matrix instantaneously conforming to the local change in interstitial melt pressure.
The model framework is applicable for relatively small mush melt extraction amounts [strains less than a few percent (Biot 1973 )] with a linear relation between ef fecti ve stress and strain. In other words, strain is assumed to be fully recoverable. Fur ther more, it is appropriate for timescales shorter than would be expected by crystal deformation due to dislocation or dif fusion creep. An y poroviscous or poroplastic effects due to rearrangements of crystals on short timescales are neglected, under the assumption that the simulated porosity changes are small enough to allow this approximation (Marchal et al. 2013 ), even if they were to occur at timescales that are comparable to poroelastic deformation (Marchal et al. 2013 ). Including such processes would require adjusting the stress balance be yond Darc y's law.
The relative importance of poroelastic, poroviscoselastic, or poroplastic behaviour will likely depend on how coherent the mush framework is. Sintering of crystals, as suggested by Holness ( 2018 ), ma y fa vour a more coherent mush skeleton, which would favour poroelastic melt extraction, whereas less coherent frameworks may allow for poroviscous or poroplastic deformation through crystal rearrangement. These are all possibilities that could be investigated in subsequent work, assuming there are sufficient constraints on constitutive relations. We note neither foliation nor lineation textures are present in the enclaves at Cord ón Caulle, which indicates a lack of evidence for significant crystal rearrangement due to settling or compaction (W inslo w et al. 2022 ). We also note microtextures that indicate recrystallization or diffusion/dislocation creep (viscous deformation of the crystals) are often not observed in magmatic systems (Holness 2018 ).
The sequence of events for the model is as follows (Fig. 3 ). An eruption occurs over some time and removes a prescribed amount of material from the rhyolite layer. The rhyolite is assumed to have a sill-like geometry and surrounded by elastically deformable rock. Therefore, a change in volume of the sill, either due to magma erupti ve withdraw al, expansion/contraction of the underlying mush, and/or melt addition from the mush, results in a change in pressure of the magma within the sill. Because the sill is assumed to be in direct contact with the underlying mush the pressure of the interstitial melt at the top of the mush must be equal to that of the overlying rhyolite. Thus, the rhyolite sill acts as a time-dependent pressure boundary condition for the underlying mush. A decrease in this pressure results in pressure gradients that induce upward flow of interstitial melt within the mush and from the mush into the overlying rhyolite.
Our model calculates how the melt pressure within the mush evolves over time and how melt flows from mush into the overlying rhyolite, as well as the resultant change in pressure of the rhyolite as it is being recharged from the mush.
Although deformation patterns at Cord ón Caulle allow for the possibility of multiple and somewhat interconnected magma bodies (Jay et al. 2014 ), co-and post-eruptive deformation can be reconciled with a sill-like magma body (Delgado et al. 2016 ;Delgado 2021 ). Based on these results, we approximate the rhyolite as an idealized rectangular sill of uniform thickness and laterally invariant magma pressure. We further assume the underlying crystal mush is contiguous with the rhyolite sill, of the same lateral dimensions, and with laterally invariant properties, including the pressure of the interstitial melt.
Our assumptions allow us to approximate the mush-rhyolite system as 1-D. Although this is a simplification, there are insufficient observational constraints to construct a geometrically more complex model. Fur ther more, the 1-D approximation enables us to focus on the key elements for our assessment of whether the post-eruptive inflation at Cord ón Caulle could be due to melt se gre gation from an underlying crystal mush. As already discussed, the assumption of no magma recharge into the mush, that is a no-flow boundary at the base of the mush, is also a simplification which could be eschewed in future investigations.
Modelling the poroelastic mush encompasses the elastic deformation of the crystal matrix and the flow of interstitial melt, both in response to a change in pressure of the overlying crystal-poor rhyolite. Interstitial melt flow results in the diffusion of pore pressure, whereas the crystal matrix responds elastically to any pressure changes. The diffusion of pore pressure constitutes the porous part of the poroelastic behaviour of the crystal-rich magma mush. The elastic part is encapsulated by the uniaxial specific storage coefficient, S s , the volume of melt released from a volume of mush per unit change in pressure, assuming a constant vertical stress and negligible horizontal strain (see Appendices A and B for detailed deri v ations). It encapsulates the compressibility of both the crystal matrix and the interstitial fluid, which is mostly rhyolitic melt, but based on evidence from the mafic enclaves also contains small bubbles of e xsolv ed volatiles (W inslo w et al. 2022 ). We therefore use the melt-bubble mixture properties throughout our analysis.
Melt flow within the crystal-rich mush follows Darcy's law
where q is the volumetric flux of interstitial melt, p is the excess (non-hydrostatic) pressure of the interstitial melt, η is the melt viscosity, k is the permeability of the mush, c = k/η S s is the ef fecti ve dif fusi vity, ζ is the volumetric change in the increment of interstitial melt content and m is the corresponding mass change (Biot 1956 ;Biot & Willis 1957 ;Wang 2000 ). Fur ther more, ρ 0 is the interstitial melt density in the reference state (Rice & Cleary 1976 ), which for our purposes will be the initial state of the model. Pore pressure and interstitial melt content are related through
where the ef fecti ve dif fusi vity, c , is assumed constant.
There are multiple ways to set up eqs ( 1 ) and ( 2) to model deformation in a magma mush. For instance, one could pose the equations in terms of mush pore pressure. An advantage of using pore pressure is the intuitive approach one can take to deriving boundary conditions. Ho wever , when pore pressure is used, eq. ( 2) becomes inhomogeneous due to changes in the mean stress with time [see Roeloffs ( 1996 ) and discussion surrounding eq. 4.63 in Wang ( 2000 )]. This adds an extra term that must be solved at every time step in the model. Using m relegates this behaviour to the boundary condition and makes the partial differential equation homogeneous. This methodology is also consistent with recent work that uses a linear poroelastic model to model deformation and magma transport within a magma mush (Liao et al. 2018 ). We therefore pose the problem in terms of m . In the following sections, we will derive the boundary conditions in terms of pressure and how to cast those conditions in terms of m .
The sill-like geometry that allows for the 1-D geometry of our model equations results in a model where pressure changes within the rhyolite are a linear function of the assumed reservoir volume. Moreover, surface uplift is calculated using a linear transfer function that is dependent on pressure. Magma withdrawal or inflow from the underlying mush will result in a mass change, M , of this rhyolite and a subsequent volume change, V , due to the displacement of the surrounding/overlying country rock, which we assume to deform elastically. Approximating the opening/closing displacement as laterally uniform for the entire sill, we reduce V to a change in thickness defined as u cu h = V / A , where A is the areal extent of the sill, u c is the displacement of the upper rock-rhyolite interface and u h is the vertical displacement of the top of the mush. Using standard formulations for displacement and pressure changes (Okada 1992 ), we define two transfer functions
P is the reduced pressure of the rhyolite, defined as P = P( t)-P ref .
P ref is the reference pressure of the system prior to eruption and assumed lithostatic and P is the total pressure. We assume that P( t = 0) = P ref .
The u transfer function allows us to calculate displacements of the host rock-rhyolite interface, u c , due to a change in reduced pressure, P , in the rhyolite. This displacement is then propagated to the surface according to w , which gives the ratio of surface displacement to u c . Fur ther more, w is the maximum vertical surface displacement (Fig. 3 ) to be compared to the LOS proxy data. We note w will al wa ys be less than u c (Okada 1992 ), since the vertical displacement from u c will translate to both vertical displacement at the surface as well as some lateral deformation.
To calculate the transfer functions, we use a MATLAB R wrapper from Toda et al. ( 2011 ) of the original FORTRAN code that solves the equations presented in Okada ( 1992 ). The sill model returns the displacement and its deri v ati ve at a point in space given a sill geometry (cross sectional area, depth, aspect ratio and dip), a shear modulus of the host rock, a Poisson's ratio, and a sill opening amount. The model then calculates the stresses and displacements at the point specified. We use this sill model to iterati vel y calculate the pressures that produces a specific change in sill thickness and associated maximum surface displacement. We then determine the linear relation between pressure and sill opening, which we define as u . The term w , in turn, is the linear relation between the sill opening and the maximum surface displacement. The value of u multiplied by the host rock's shear modulus is a constant, and w is constant for a given sill geometry.
The change in mass of the rhy olite lay er, M , leads to a Dirichlet pressure boundary condition on the mush. That change can be expressed in two ways: through the rhyolite fluxes in and out of the layer, and by the changes in volume and density of the melt. In terms of fluxes,
with M mush being the rhyolite mass leaked from the mush and M E being that erupted from the rhyolite layer. This difference in mass fluxes is then expressed in terms of Q , the rate of mass erupted (e.g. Fig. 4 ) and q h , the Darcy flux at the top of the mush from eq. ( 1). Details of its deri v ation can be found in Appendix B2.2.
Given that M = ( ρV ), the time rate of change in M can also be expressed as
where following Rice & Cleary ( 1976 ), we have linearized ρ and V .
Here, h r is the initial thickness of the rhyolite layer. Throughout, the subscript 'r' denotes properties of the rhyolite layer. Fur ther more,
is the compressibility of the rhyolite melt plus v olatile bubb les mixture, which we assume constant (e.g. Rice & Cleary 1976 ;Anderson & Segall 2011 ). We define the inverse bulk modulus of the interstitial melt-bubble mixture, 1/ K r , by assuming an ideal gas phase for the supercritical fluid mixture inside bubbles
Here, φ g is the volume fraction of bubbles in the mixture. One can use an equation of state other than an ideal gas, but considering the uncertainties in the model parameters which are discussed in Section 5 and Appendix C1, we assume ideal behaviour. A final form of eq. ( 5) that includes an expression for d u c /d t from eq. ( 3) is
Equating the two expressions for the time rate of change in mass of the crystal-poor rhyolite (eqs 4 and 8 ) gives
The time rate of change for the displacement of the top of the mush,
, implicitly depends on P . An expression for it is derived in Appendix B2.3, and is
Combining eq. ( 10) with eq. ( 9) gives
where
is the scaling between change in length and change in pressure of the crystal-poor rhyolite magma, β eff is the effective compressibility of the crystal-poor rhyolite magma, and
is the volumetric magma eruption rate. Each term in β eff (eq. 12 ) corresponds with a compressibility for different components of the crystal-poor rhyolite magma. The first term, 1/ K r , is the compressibility of the mixture of the rhyolite melt and v olatile bubb les, u / h r represents how compressib le the host rock is, and h/ h r ˜ K v defines the strength of the top of the mush. Interestingly, this last term modulates the undrained compressibility of the mush by the ratio of the mush size to the rhyolite size. Relati vel y larger mushes will impact the ef fecti ve compressibility through this term compared to larger rhy olite lay ers. The term then scales this ef fecti ve compressibility to a change in length given a change in pressure. The three terms 1/ K r , u / h r and 1 / ˜ K v are each a component of the rhy olite lay er's strain due to compressibility, so that the change in length of the rhyolite is the sum
We assume that our simulations start from a static initial state with no flow of interstitial melt within the mush. Thus, at t = 0 we assume m ( z ) = 0 and ∂ p / ∂ t = 0. This is likely a simplification, given that prior to the eruption there was inflationary deformation, albeit at a significantly lower rate than deformation during and immediately after the eruption (Delgado 2021 ). Using a non-stationary initial state would require a different and more complicated set of assumptions. We therefore choose the simpler approach of assuming a static initial state. Future w ork, ho wever , could attempt to assess the effect of a non-stationary initial state.
Eruption rates for the Plinian to subplinian phase of the eruption have been estimated on the basis of plume height (Fig. 4 ) to range from approximately 3 × 10 6 to 3 × 10 7 kg s -1 initially and subsequently decreased to 10 5 to 10 6 kg s -1 (Bonadonna et al. 2015b ).
The subsequent ef fusi v e phase, which included vigorous ash v enting (Schipper et al. 2013 ;Crozier et al. 2022 ) and formation of a lava flow, had initial effusion rates of approximately 3 × 10 5 kg s -1 that decreased over time to approximately 3 × 10 4 kg s -1 . These rates are estimated from satellite thermal data (Coppola et al. 2017 ) and from volume changes of the lava flow that were estimated from satellite images taken at different times throughout the eruption (Bertin et al. 2012(Bertin et al. , 2015 ) ).
A function for the mass eruption rate Q ( t ) is required for the calculation of P (eq. 11 ). We obtain Q ( t ) from a smoothing cubic spline interpolation of the aforementioned discharge rate estimates shown in Fig. 4 . Calculation of the smoothing spline was done using the MATLAB R function csaps , which minimizes the error between Q ( t ) and the discharge rate estimates under a smoothness constraint, θ. The data w ere w eighted to obtain a monotonically decreasing function of Q with respect to time. The value of Q ( t = 0) was chosen such that the time integral of Q falls within a reasonable range of estimates for the total erupted mass during the first few days of the eruption (Bonadonna et al. 2015a ;Pistolesi et al. 2015 ), which are shown in Fig. 4
As already discussed, at the bottom of the crystal mush, we assume a no-flow (Neumann) boundary condition
In reality a shallow crustal magmatic system may be contiguous at depth with a more e xtensiv e crustal magma transport and storage system. If so, one question would be to what extent permeability and other proper ties var y with depth. Moreover, it is hard to imagine that the mush is not an open system in the sense that there would not be a deeper magma supply, at least on longer timescales (e.g. Cashman et al. 2017 ;Hildreth 2021 ). Because our objective is to assess to what extent the obser ved inflationar y defor mation since the eruption can be explained by melt percolation from a crystal mush, as opposed to recharge, we assume a no-flow boundary at the bottom of the mush. In other words, some extent of magma recharge from below can al wa ys be invoked to explain volcano inflation, here we wish to assess the viability of the closed-system response on the decadal timescale of a potential mush beneath Cord ón Caulle. It should be noted that the no-flow assumption is also consistent with the aforementioned static initial state approximation, allowing us to isolate the response of the magma mush to the eruption. At the top of the crystal mush ( z = h ), that is at the interface between the mush and overlying crystal-poor rhyolite melt, the pressure of the interstitial rhyolite melt within the mush and within the crystal-poor rhyolite must be continuous. We thus impose a time-dependent Dirichlet boundary condition that represents the reduced pressure, P , of the crystal-poor rhyolite melt in contact with the top of the mush. The value of P is calculated using eq. ( 11), with the added complication that this pressure boundary condition must be expressed in terms of m , for use with eq. ( 2). Therefore, P must be expressed in terms of m . This relation is given by (see Appendix B2.3)
where m h = m ( z = h ). Combined with eq. ( 11), eq. ( 15) is the boundary condition for the top of the mush, and the entire problem to be solved is given by eqs ( 2 ), ( 11), ( 14) and ( 15 ).
To predict surface displacement, w , we incorporate the transfer functions (see Section 4.3.1 for a more detailed discussion of how the transfer functions are calculated) as calculated by Okada ( 1992 ) in eq. ( 3) and write
The term w is the predicted maximum vertical surface displacement, which occurs directly over the centre of the sill. The LOS displacement may be different from this value, depending on the angle of incidence, the direction of deformation, and where on the volcano it is measured. Our LOS proxy data (Fig. 1 b) has some variability in incidence angle and location (Delgado 2021 ). We also do not know the direction (compared to vertical) of maximum deformation for these LOS proxies. This would require re-analysing the descending and ascending interferograms that led to the LOS points to retrieve the vertical and horizontal components of the LOS displacement.
Alternati vel y, one could assume a pseudo-vertical displacement by di viding the LOS b y the cosine of the incidence angle (e.g. Delgado et al. 2016 ), although this assumes all displacement to be in the vertical direction, which would likely overestimate the actual vertical displacement. For simplicity, we therefore compare w to the LOS displacement. Also, our objective is to assess the feasibility of melt se gre gation from a mush causing uplift. Comparing w to the LOS displacement directly allows us to test this feasibility.
Eq. ( 2), together with boundary conditions given by eqs ( 11) and ( 14), as well as eq. ( 15), are solved by the Method of Lines (e.g. Fletcher 1998 ) using the MATLAB R solver ode15s . The main result of these simulations is the prediction of m ( t ), which allows the estimation of the surface uplift, w , using the relations from eqs ( 15) and ( 16). The values of w ( t ) can be compared against the LOS displacements shown in Fig. 1 (c).
The model defined in the previous section, comprised of eqs ( 2), ( 11), ( 14), ( 15) and ( 16), depends on six poroelastic parameters: K s , G , φ 0 , α, k and S s , plus melt viscosity η. The remaining poroelastic parameters are dependent on these parameters as summarized in Table 1 . Parameters that do not directly affect the mechanics of the poroelastic mush, but ne vertheless af fect the simulation results, are summarized in the lower part of Table 1 and include the mush thickness, h , the thickness of the rhyolite, h r , the mush depth, d , the areal extent, A , of the Cord ón Caulle magmatic system and the shear modulus of the host rock, G c . Most of these parameters can be constrained based on prior information. Only three parameters, S s , h and h r , are truly unknown and will be estimated through model simulation.
To achieve the parameter estimation, we solve the aforementioned equations in order to predict vertical surface displacements that we match to the LOS displacements (Fig. 1 c) as closely as possible. As discussed pre viousl y, this approach neglects the fact that the LOS displacements are not exactly equal to vertical surface uplift. Both the poroelastic response of the mush and the eruption rate Q jointly result in surface deformation. Because the eruption lasted for almost one year, the response of the mush to eruptive magma withdrawal from the overlying rhyolite already begins during the eruption. It is therefore difficult to limit our simulations solely to the time after the er uption. Fur ther more, parameters that do not govern the response of the mush nevertheless affect model predictions, in particular the magnitude of predicted surface deformation, both during deflation and inflation. Therefore, we simulate the entire time period, starting from the beginning of the eruption until early 2020.
To quantify the match between simulated surface deformation and the LOS displacements we define a model misfit, χ , which is the root mean square error between the simulated surface deformation and the LOS displacements. It is given by
Here, LOS j corresponds with the j th LOS displacement value from Delgado ( 2021 ) and w j is the corresponding simulated (predicted) maximum vertical displacement based on the assumed sill geometry (discussed in Section C1.8 and shown in Fig. 1 b). The weights W j Table 1. Model parameters related to the poroelastic mechanics of the crystal mush and their interdependencies as well as other parameters that define the geometry of the reservoir and the rock properties. See Appendix C1 for explanations of the parameter ranges and Section 4.3.1 for an explanation of how the transfer functions are calculated.
Units Definition Value Poroelastic parameters α (-) Biot-Willis coefficient 0.5 -0.9 φ 0 (-) Initial mush porosity 0.36 K s (Pa) Bulk modulus of the crystals matrix 87.6 × 10 9 G (Pa) Shear modulus of mush 11.2 × 10 9 η (Pa s) Melt viscosity 3300 k (m 2 ) Mush permeability 1.7 × 10 -9 -6.7 × 10 -9 a (m) Grain size 0.001-0.002 S s
(Pa -1 ) Uniaxial specific storage coefficient
Rhy olite-bubb le mixture bulk modulus
Vertical drained bulk modulus of mush
Pressure dif fusi vity of the mush k /( ηS s ) Other parameters h (m) Initial mush layer thickness Estimated h r (m) Initial rhyolite layer thickness Estimated A (m 2 ) Magma chamber areal footprint 60 × 10 6 m 2 h r A (km 3 ) Rhy olite v olume > 1.65 km 3 d (m) Magma chamber depth 6000 m G c (Pa) Host rock shear modulus 1 -35 × 10 9 Pa u (m Pa -1 ) Transfer function for pressure to rock displacement 4621.7/ G c w (-) Transfer function for rock displacement to surface displacement 0.3324 (m Pa -1 ) Scaled compressibility of the magma chamber h
have a value of one, except W 1 and W 2 , which are zero because of the large rate of change in LOS during the first few hours of the eruption as well as the large uncertainty in eruption rates during that time (Fig. 4 ). Our weighting of the LOS data ensures that low values of χ meet our objective of quantifying the misfit to the inflationary part of the LOS data.
Although we are searching for the best fit to the LOS displacements, we wish to reiterate the inherent imperfections in this data. The LOS proxies are a compilation of varying orbital angles and a moving locus of maximum defor mation, par ticularly during the eruption, which have been spliced together to provide a continuous signal. As such, the exact misfit value should be viewed with these caveats in mind. We also note that our goal is not to match the uplift data any better than was already done using elastic and viscoelastic models (Delgado et al. 2016(Delgado et al. , 2018 ) ). We only wish to illustrate what is necessary for a mush to produce all of the measured uplift.
Because the dimensionality of parameter space makes a grid search prohibitive, we conducted an exploratory random sampling (e.g. Kochenderfer & Wheeler 2019 ) of all parameters within some predefined range. The exact parameters values and ranges as well as dependencies are summarized in Table 1 . Explanations of how we came to these ranges and values are given in Appendix C1. Random sampling herein refers to using MATLAB R 's pseudo-random number generator to produce a set of values for each parameter within the range we wish to cover. The purpose of this random sampling was to explore parameter space with respect to relationships between parameters. Our results of this random sampling are presented in Appendix C1.12.
We found that for any combination of parameters there is a unique value of both S s and η/ k for which the misfit χ is the smallest. This is denoted as χ opt . For some combinations of parameters the smallest value of χ found within S s -η/ k space corresponds to the minimum in χ that can be achieved for any combination of parameters. We denote this value as χ min . We emphasize that it is not a unique minimum, because more than one combination of parameters can produce this minimum value of χ . We seek parameter combination that result in χ min .
The parameter combination A ( α = 0.7, K s = 88 GPa, G = 15 GPa and φ 0 = 0.42) gives the minimum misfit value of χ opt = χ min = 3.75 cm (Fig. 5 a). In contrast, the combination F, comprising a higher value of α, a lower value of φ 0 and a lower value of G , gives a corresponding smallest misfit value of χ opt > χ min (Fig. 5 b). The modelled solution corresponding to A provides a distinctly better fit to the LOS time-series than that corresponding to F (Fig. 5 c).
Four other simulations, labelled as B, C, D and E in Fig. 5 (d), illustrate how changes in parameters affect the misfit. Parameters can affect the rate at which the poroelastic response occurs, the magnitude of the response, or both, where the response magnitude is the final LOS displacement amount. The rate at which the mush responds to pressure changes in the rhyolite is dependent on the characteristic pore-pressure diffusion time of the mush, τ = h 2 S s η/ k . Because S s and h also affect the response magnitude in simulated uplift, a change in η/ k will influence the rate of uplift. For instance, higher permeability enables more rapid melt segregation, facilitating an accelerated poroelastic recovery. The final amount of recovery remains unchanged, leading to the same final LOS displacement value as a lower permeability case. Ho wever , the extent of deflation will vary due to the faster recovery, which corresponds to increased melt se gre gation during the eruption and leads to less deflation. Varying any of the parameters A , h r , G c , u and w will change the magnitude of simulated uplift only, shifting the entire curve vertically. Here it should be kept in mind that the simulated uplift is also dependent on Q , the mass eruption rate, which w e ha ve not varied betw een simulations from what is shown in Fig. 4 .
The values of χ opt , which are obtained from a random sampling of all parameters with minimization over S s and η/ k , can be shown to ha ve w ell defined functional dependencies (see Appendix C1.12). These functional dependencies are expressed in terms of three functions 1 , 2 and 3 that encompass the three unknown parameters S s , h and h r (C2). It is therefore possible to leverage these functional relations to find which parameter combinations yield the minimum misfit of χ min = 3.75 cm. Specifically, for a given set of a priori parameters that are either constrained to a specific value or to within a defined range (see Table 1 ), we can find those combinations of S s , h and h r that have χ as the smallest misfit in S s -η/ k space.
The best fits of model simulations to the LOS displacements, that is solutions corresponding to χ min , provide estimates for the specific storage of the mush, S s , the thickness of the mush, h and the thickness of the overlying crystal-poor rhyolite, h r . As already discussed in the previous section, these estimates are based on prior constraints of other parameters, as summarized in Appendix C1 and in Table 1 . Fur ther more, it should be kept in mind that the LOS displacements provide a proxy for deflationary and inflationary deformation at Cord ón Caulle, whereas the simulated time-series represents vertical displacements that are from a model that is based upon a number of simplifying assumptions. Thus, the results presented here should solel y be vie wed within the context of testing the the feasibility that inflationar y defor mation after the er uption could be due to melt se gre gation from a crystal mush. In other words, our parameter estimates provide an assessment of the extent to which the conditions required for the crystal mush hypothesis to hold are feasible.
Throughout, we assume a homogeneous mush with constant values for all parameters. In nature this seems unlikely. For example, a magma mush beneath Cord ón Caulle could be stratified such that porosity, which affects k , or volatile bubbles within the interstitial melt, which affects S s , could vary with depth (e.g. Bachmann & Bergantz 2004 ;Hildreth & Wilson 2007 ;Cashman et al. 2017 ;Sparks et al. 2019 ). If that were the case, then our model would be more representative of the shallower parts of a vertically more extensive mush system (e.g. Hildreth 2004 ;Hildreth & Wilson 2007 ;Cashman et al. 2017 ;Hildreth 2021 ).
We present our results in Fig. 6 as a function of mush permeability, k , which is treated as if it were an independent parameter, even though we have constraints on k that are based on the Cord ón Caulle enclaves (see discussion in Appendix C1). The range is based on values that are consistent with permeabilities for partially crystallized rocks and other porous materials at similar porosities and grain sizes as the enclaves (Hersum et al. 2005 ). Ho wever , because it is well known that field-scale permeabilities can be considerably larger than permeabilities measured at the sample-scale (e.g. Clauser 1992 ;Schulze-Makuch et al. 1999 ), for example due to the presence of highly permeable pathwa ys, w e ha v e e xtended the permeability range considered by an order of magnitude.
We estimate the specific storage of the mush, S s , to fall within the range of about 10 -10 to 10 -9 Pa -1 (Fig. 6 a). The specific storage represents the amount of interstitial rhyolite released from the mush per unit volume of mush. It depends on the elastic deformability of the crystal matrix of the mush, as well as the compressibility of the pore fluid, in this case rhyolite melt with v olatile bubb les. The value of S s therefore is strongly dependent on the volume fraction of bubbles within the interstitial melt, as indicated by eqs ( 7) and ( C3 ). It is commonly assumed that mush systems associated with evolved magmas contain exsolved volatiles (e.g. Wallace 2001 ;Shinohara 2008 ), due to crystallization-driven exsolution ( second boiling ) and/or the degassing of deeper recharge of magma. As discussed by Parmigiani et al. ( 2016 ), despite their b uoyancy, b ubbles can remain trapped within a mush due to capillary and viscous forces at values of φ g ≤ 0.15 (see discussion in Appendix C1). We therefore consider solutions of S s ≤ 4.3 × 10 -10 Pa -1 as viable, as indicated in Fig. 6 (a).
Within the observational constraints of enclave permeability and v esicularity, the observ ed uplift can be e xplained by melt se gre gation from a mush. Where these conditions are all met is indicated by the circle and arrow in Fig. 6 . The permeability here is k ≈ 6.7 × 10 -9 m 2 , and the corresponding specific storage is S s ≈ 4 × 10 -10 Pa -1 , which corresponds to φ g ≈ 0.15. Permeabilities of k < 6.7 × 10 -9 m 2 , which fall within the range of sample-scale estimates obtained from measured enclave porosities and crystal sizes. Because they require gas fractions of φ g > 0.15 we consider them as unlikely. In contrast, the possibility that the effective permeability of the mush at the field-scale is somewhat larger than the samplescale estimates would allow a larger volume of mush with values of S s corresponding to gas fractions as low as a few percent.
The corresponding estimates for the crystal-poor rhyolite thickness (Fig. 6 b) and underlying crystal mush (Fig. 6 c) are h r ≈ 260 m (volume of V r ≈ 15.6 km 3 ) and h ≈ 1 km (volume of V ≈ 60 km 3 ), respecti vel y. To put these values in perspective, the predicted crystal-poor rhyolite volume is ten-fold of what erupted in 2011-2012, implying that more than 14 km 3 of crystal-poor rhyolite could still remain beneath Cord ón Caulle. The rhyolite predicted to be contained within the mush is slightly more than that. Approximately 25 per cent of the erupted volume would have been recharged by melt se gre gation from the mush, which equates to a within the sample-scale permeability estimated from the enclaves and with φ g no greater than the value at which bubbles are expected to be trapped within a crystal mush (Parmigiani et al. 2016 ). The maximum porosity change is expected to be 1 per cent here, which is well within the range of poroelastic behaviour (e.g. Karig & Hou 1992 ). Lower permeabilities than this value that meets all r equir ements require higher gas fractions to produce the proper rhyolite compressibility than are feasible. This is not the case for higher permeabilities, if one considers them viable based on the consideration of sample-versus field-scale permeabilities. All results shown are for a host rock shear modulus of G c = 12 GPa, which equates to an eruptive pressure change of 32 MPa, and an areal footprint of A = 60 km 2 . We allow α to vary with dark blue representing α = 0.7, solid light blue as α = 0.9, and dashed light blue as α = 0.5. v olume of appro ximately 0.4 km 3 . Thus, appro ximately 21 km 3 of rhyolite would remain in storage. The aforementioned estimates are somewhat larger or smaller if one assumes a different value for the Biot-Willis coefficient (Fig. 6 ).
We also estimate the maximum change in mush porosity, φ, in Fig. 6 (d), which is the absolute maximum change in porosity over the course of the simulation (see Appendix B3). A maximum increase of 1 percent is anticipated under the aforementioned constraints. This falls within the range of porosity changes observed during recoverable deformation of unconsolidated granular materials with similar porosities as the enclaves (Karig & Hou 1992 ). Porosity changes are inversely correlated with permeability. At permeabilities less than ≈2 × 10 -9 m 2 , predicted porosity changes are likely too large for the poroelastic assumption to hold. Thus, both the required vesicularity and the predicted porosity changes make solutions corresponding to permeabilities less than 6.7 × 10 -9 m 2 unlikely. Instead, higher permeability solutions are preferred, if field-scale permeabilities are applicable.
We adopt an areal extent of the Cord ón Caulle system of 60 km 2 . Both smaller and larger areal footprints have been proposed (e.g. Delgado et al. 2016 ;Novoa et al. 2022 ). Smaller areas would result in minimum misfit combinations for sample-scale permeabilities that violate the maximum φ g constraint of 0.15. If permeabilities are indeed larger that 6.7 × 10 -9 m 2 due to field scale effects, then smaller areas could still have viable solutions. Larger areal extents push the results further into the allowable regime, meaning smaller permeabilities ( < 6.7 × 10 -9 m 2 ) would then be feasible. We chose the current geometry on the basis that the magma storage system would not extend outside of the Cord ón Caulle graben.
The simulations were for a shear modulus of the host rock of 12 GPa, which is within the range of values used in previous studies of deformation at Cord ón Caulle (Delgado et al. 2019 ;Novoa et al. 2022 ). The simulated deflation during the eruption consequently corresponds to a decrease in pressure of the crystal-poor rhyolite of 32 MPa, which is similar to the overpressures estimated by Novoa et al. ( 2022 ). The corresponding pressure increase associated with the inflationary deformation after the eruption is 7 MPa. For different shear moduli of the host rock the simulated pressure changes would also be dif ferent. Howe ver, for smaller shear moduli, the estimated values for S s would require values of φ g that violate the aforementioned upper limit of 0.15.
Our simulations result in a smoothly decaying rate of inflation after the eruption. The LOS displacements (Fig. 1 c) show some low amplitude higher order temporal variability superimposed on this long-term trend, which our model cannot explain. Assuming that these variations are indeed magmatic, then some additional processes would need to be invoked to explain these small variations on top of the broader inflationary trend. They could include, for example, transient changes in the hydraulic properties of the magmatic system (e.g. Le M ével et al. 2021 ), temporal variations in volatile content (Spang et al. 2022 ), or several diking episodes. For the latter there is however no seismic evidence (Delgado et al. 2016(Delgado et al. , 2018 ) ).
Are the estimated volumes of the mush and overlying rhyolite reasonable, and do they support the potential existence of a crystal mush beneath Cord ón Caulle? If taken at face value, meaning that our estimates encompass the entire magma storage system beneath Cord ón Caulle, they would be consistent with what has been inferred by Singer et al. ( 2008 ). Based on their synthesis, Cord ón Caulle has erupted e xclusiv ely rhyolitic and rhyodacitic magmas since 32 ka, with half of the rhyolitic activity during the last 16.5 ka. In particular, the 1921In particular, the -1922In particular, the , the 1960In particular, the , and the 2011In particular, the -2012 rhyolites represent some of the most evolved magmas erupted in the entire Puyehue-Cord ón Caulle complex, although it should be noted that the extent to which Puyehue and Cord ón Caulle represent the same system remains unclear (Singer et al. 2008 ). Prior to the past 32 000 yr, that is between 173 and 32 ka, magmas were mainly basaltic and andesitic in composition, with a gap of > 9 wt per cent in SiO 2 separating the older lavas from the later silicic ones. Jicha et al. ( 2007 ) found that the rhyolites produced by the entire complex may have fractionated from 85 per cent crystallization of a basaltic to a basaltic andesitic parent magma. Consequently, about 130 km 3 of basaltic magma would have been required to produce 14.5 km 3 of rhyolitic magma through fractional crystallization (Singer et al. 2008 ). We find a mush volume of 60 km 3 , if we allow for sample-scale permeabilities only, which is about half of the total volume of basalt predicted by Singer et al. ( 2008 ). Ho wever , the possibility of field-scale permeabilities that are larger than permeabilities at the sample scale, or vertical stratification in mush properties, would result in volume estimates of 130 km 3 or more.
If one takes a broader view of crystal mush systems which erupt crystal-poor silicic magmas it is feasible, albeit speculative, that our anal ysis onl y captures the upper parts of a verticall y hetero genous or stratified system that has been evolving as a consequence of deeper mafic magma supply and gravitationally driven segregation processes (Bachmann & Bergantz 2004 ;Hildreth 2004 ;Cashman et al. 2017 ). For example, one could speculate that the Cord ón Caulle system at present is similar to the Long Valley system during its early Glass Mountain stage (Hildreth 2004(Hildreth , 2021 ) ). This inv ariabl y leads to the question of long-term thermal viability and evolution of the Cord ón Caulle system. Hildreth ( 2021 ) considers that the longevity of the melt percolation process is favoured by several kilometres of permeable mush with hot mafic recharge from below. At the same time he leaves open w hether rhy olite melt bodies that form atop a mush are thermall y buf fered for long time periods or whether they are intermittent features. Our model neither provides insight into the longer-term evolution of such a system, nor does it include the role of deeper magma supply. This is by design. The no-flo w lo wer boundary condition for our model was chosen to provide the simplest model possible that can represent the shortterm poroelastic response to eruptive magma withdrawal from the crystal-poor rhyolite.
A finite mush is also required by the poroelastic model because the rate of LOS uplift exponentially decays over time. Within the constraints of a poroelastic model and depth-invariant properties, this can only be achieved with a finite mush where the pore-pressure diffusion timescale is similar to the exponential decay timescale of the LOS uplift rate. In other words, the change in pore pressure that dif fuses downw ard through the mush reaches the bottom noflow boundary, and thereby causes the decay in the rate of observed uplift. For a mush that is suf ficientl y thick so that the pressure change does not reach the bottom boundary, the rate of uplift remains nearly constant throughout the simulation period, unless there is a decrease in permeability or increase in specific storage with depth. Here we assume mush properties that are spatiall y inv ariant, reco gnizing that any changes in properties with respect to depth would occur over lengths that scale with the estimated mush thicknesses shown in Fig. 6 (c).
Our results are for predicted changes in porosity of the order of 1 per cent (F ig. 6 d), w hich is consistent with the observed poroelastic behaviour of unconsolidated granular materials (e.g. Karig & Hou 1992 ). As already discussed in Section 3, the potential rele v ance of poroviscoselastic or poroplastic beha viour ma y depend on the degree of sintering of crystals (Holness 2018 ), which provides impetus for further investigation of samples for evidence of such deformation. At present such evidence seems lacking (W inslo w et al. 2022 ;Holness 2018 ).
The existence of the mafic enclaves with interstitial rhyolite points tow ards a primiti ve mush (W inslo w et al. 2022 ). Much of the 85 per cent of chemical fractionation from crystallization of a basaltic to a basaltic andesitic magma envisioned by Jicha et al. ( 2007 ) is evident in the enclaves (W inslo w et al. 2022 ). Since this primitive mush is presumably producing the rhyolite, little fractionation from deeper primitive melts is necessary. Thus, substantially deeper melt storage is not required to explain the recent events at Cord ón Caulle. Ho wever , we do not discount recharge as an important component in this magmatic system, because the long-term viability of the system will of course depend on this (e.g. Gelman et al. 2013 ;Annen et al. 2015 ). In ter ms of ther mal viability and long-term evolution of magma-mush systems (Gelman et al. 2013 ), we point to the existence of a large geothermal system at Cord ón Caulle (Sep úlveda et al. 2005 ) for which Singer et al. ( 2008 ) proposed latent heat release from crystallization as a partial driver during the past 16.5 ka.
Although we argue that the observed uplift at Cord ón Caulle is a response to melt se gre gation from a magma mush, no direct geophysical observation of this process has been recorded (although the existence of the enclaves lend credence to this hypothesis). Our model is one of a few physical models that could explain the uplift, including models that invoke deeper seated recharge (Delgado et al. 2016 ) or viscoelastic relaxation of the host rock (Segall 2016 ). Like the model we have presented, viscoelastic relaxation also may not require recharge to produce uplift in specific circumstances (Segall 2016 ).
To differentiate between these possible mechanisms, future geophysical surv e ys could be completed, as each mechanism ma y ha ve distinct signatures. For example, imaging the magma reservoir may be possible using gravity surv e ys (e.g. Flinders et al. 2013 ;Harpp & Geist 2018 ), magnetotelluric measurements (e.g. Samrock et al. 2021 ) and seismic tomography (e.g. Lees 2007 ), or combinations of these techniques (Pritchard & Gregg 2016 ), assuming a high enough resolution could be obtained. Results might show whether a mush exists beneath Cord ón Caulle, possibly constrain its size, and also provide more information about the architecture of the greater Puyehue-Cord ón Caulle system. Although these future measurements may help constrain deformation models, they cannot necessarily distinguish between one mechanism over another. For instance, the presence of a mush does not necessitate that poroelastic ef fects completel y contribute to surface uplift. Ho wever , the lack of a mush would rule out this mechanism.
Temporal observations may provide more insight into underlying processes ongoing within the reservoir. Temporal changes in gravity may help determine causes of ongoing uplift (Battaglia et al. 2003 ;Segall 2010 ). Recharge from deep should have a stronger positive gravity variation compared to either mush melt migration or viscoeleastic relaxation. Recharge brings new mass into the reservoir, causing a relati vel y large positi ve v ariation, whereas mush melt migration only transports material a short distance, leading to a smaller gravitational change. Viscoelastic relaxation without recharge may induce an even smaller gravitational change after accounting for the deformation, as material is generally not being transported. Temporal gravity measurements may even discriminate the contribution of each mechanism to the observed deformation. Gas and thermal monitoring of the hydrothermal system also may help differentiate causes of inflation. Deep-sourced recharge may be more volatile rich and also hotter than rhyolitic melt percolating from the mush (Edmonds 2008 ).
Intriguingl y, the post-erupti ve uplift at Cord ón Caulle has been nearly aseismic (Delgado et al. 2016(Delgado et al. , 2018 ) ). Delgado et al. ( 2018 ) invoked the Kaiser effect to explain the lack of seismicity (Heimisson et al. 2015 ), which predicts that any seismicity will remain low until the stress in the magmatic reservoir is close to the maximum value of the previous loading cycle (i.e. prior to eruption). Our model cannot help e v aluate this hypothesis. One could ask though whether seismicity would be expected during poroelastic melt migration, in contrast to magma injection, which has been envisaged alternati vel y (Delgado et al. 2016(Delgado et al. , 2018 ) ). A possible alternate model for inflation at Cord ón Caulle could be based on viscoelastic relaxation of the host rock (Segall 2016 ), although it is not clear whether this would be mostly aseismic. The uplift caused by the melt se gre gation from a crystal mush w ould, ho wever , still require the Kaiser effect at least to explain the lack of seismicity within the host rocks. Future monitoring of the seismicity at Cord ón Caulle, along with deformation, may help shed led light on the ongoing processes.
The 2011-2012 eruption at Cord ón Caulle was fed from a crystalpoor rhyolite body that should underlie a large part of the Cord ón Caulle graben. Our analysis supports the possibility that posteruptive uplift at Cord ón Caulle may be due to recharge of this crystal-poor rhyolite magma reservoir, as a consequence of melt se gre gation from an underlying crystal mush, and in response to magma withdrawal from the rhy olite reserv oir during the 2011-2012 eruption. Neither deeper recharge during the time period since the eruption ended (Delgado et al. 2016(Delgado et al. , 2018 ) ) nor alternate processes that give rise to inflationary deformation, such as viscoelastic relaxation of the host rock (Segall 2016 ), can be ruled out through our anal ysis. It is, howe ver, possible to account for all of the uplift of the past decade through melt se gre gation from a magma mush.
The potential existence of such a mush is supported by the presence of mafic enclaves within the erupted rhyolite (W inslo w et al. 2022 ). Moreover, crystal mush has been postulated as the source for evolved silicic magmas (Hildreth 2021 ), such as those that have been erupting at Cord ón Caulle during the past 16.5 ka. In addition, the release of latent heat associated with the evolution of such a crystal mush would likel y suf fice to drive Cord ón Caulle's geothermal system. Thus, multiple lines of reasoning support the potential existence of a crystal mush beneath Cord ón Caulle. If this was indeed the case, then mush-derived recharge as envisaged and simulated herein would arguably be an inescapable consequence of the eruption.
Geometrically idealized 1-D modelling of this process, together with observational constraints, results in estimates of relevant parameters. In light of the many simplifications and assumptions, these estimates should be viewed as laden with uncertainty and are solely for the purpose of testing the feasibility of the mush hypothesis. Within that context, best estimates include a permeability of 6.7 × 10 -9 m 2 for the mush, based on average grain sizes and porosities of the Cord ón Caulle enclaves. Other parameter estimates are the specific storage of the mush with a value of 4 × 10 -10 Pa -1 and a mush thickness of approximately 1 km, corresponding to a volume of 60 km 3 . The overlying rhyolite would correspondingly have a thickness of 260 m or a volume of 15.6 km 3 . The corresponding increase in pressure of the crystal-poor rhyolite reservoir would be 7 MPa since the 2011-2012 eruption ended.
This work was supported by NSF-EAR 182425 to HMG, NSF-EAR 1823122 to PR and NSF-EAR 1824160 to MEP. Discussions with Kyle Anderson, Jacob Jordan and Cin-Ty Lee are appreciated. We also thank Christian Huber and Paul Segall for their insightful and v aluable re vie ws.
The data underlying this paper and the code used to produce it will be shared on request to the corresponding author. and conservation of mass,
provides two equations that govern fluid flow through the porous medium. Darcy's law states that the flux of fluid out of a porous material ( q ) is proportional to the ne gativ e pressure gradient of the fluid. The two are related by the ratio of the permeability of the matrix to the viscosity of the fluid. The two can be combined by taking the divergence of Darcy's law and eliminating ∇ • q .
Lastly, taking eq. ( B8) and substituting in ∇ 2 p leads to a diffusion equation:
with an ef fecti ve dif fusi vity,
B2 Boundaries of the mush and initial condition
Assuming no material can cross into or out of the mush sets a no flux boundary condition. Physically, this means the bottom boundary is isolated from any injection of new material, and pore melt cannot flo w out belo w. This can be done by setting the flux to zero in Darcy's law (eq. B9 ), giving
and by substituting eq. ( B8 ) for ∇p gives an expression in terms of m ,
If there is a flux into the mush at the bottom, then q can be non-zero, but that is beyond the scope of this investigation.
After an eruption, the mass of the rhy olite lay er will change with time depending on the amount entering from the mush ( M mush ) and the amount leaving the layer due to the eruption ( M E ):
One can express M mush in terms of the change in mass of the interstitial melt content of the mush m since the net change equals the amount leaked from the mush. This is written as
where h is the thickness of the mush. The ne gativ e sign comes from the convention that positive values of m indicate material entering the mush. By substituting this into eq. ( B16) and taking the time deri v ati ve, the time rate of change of M becomes
where Q is the rate of mass erupted (e.g. Fig. 4 ). Substitution from eq. ( 2) gives
where q h is the Darcy flux at the top of the mush from eq. ( 1) with m = ρ 0 ζ . This must be multiplied by the cross sectional area ( A ) to account for the 1-D flux. The change in mass can also be written according to the time dependent changes in volume and density,
which can be linearized as
The change in volume, V , is the difference between the final and initial volumes. The initial volume is simply the initial length of the rhy olite lay er multiplied by A . The final v olume depends on how much the top of the mush and the top of the rhy olite lay er have mov ed, u h and u c , respectiv ely. Increases in the top of the mush decrease the rhyolite thickness, whereas increases in the top of the rhyolite increase the thickness. In other words,
The displacement of the top of the rhyolite is related to the pressure of the rhyolite, P , by the transfer function u ≡ d u c /d P , which depends on the geometry and the depth of the magma chamber as well as the physical properties of the rock surrounding the chamber. Thus, u c = u P . Based on the definition of the rhyolite compressibility (1/ K r ),
ρ can be approximated as
Taking the time deri v ati ve of eq. ( B21 ) and substituting these previous expressions in gives
Combining eqs ( B19 ) and ( B25 ) and solving for d P /d t gives ) This defines how the pressure will change in the rhyolite layer, depending on the flux from the mush and the displacement of the top of the mush. The next steps are to find an equation for the displacement of the mush and relate the rhyolite layer pressure to m h .
Recalling the quasistatic approximation, this gives a requirement for the top of the mush. Going back to eq. ( A20 ), σ must be a purely time dependent function. This means the stress at the top of the mush is the same as the stress in the middle of the mush. Assuming stress continuity across the mush-rhyolite boundary, σ must be at all times σ = -P . (B27)
Turning to eqs ( A20 ) and ( A21 ), equations for p and ∇ • u can be derived. First, to get an equation for p , substitute σ = -P , while solving for and eliminating ∇ • u from both equations,
To solve for ∇ • u , substitute eq. ( B28 ) into the first equation of the previous sequence,
To simplify the last term, recall S s = α 2 / K v + S ,
Substituting this back in,
Integrating eq. ( B31) with respect to z over the thickness of the mush leads to
Then, taking a deri v ati ve with respect to time,
Recall that P only depends on t , so it only is af fected b y the time deri v ati ve and not the integration. Since ∂ m / ∂ t is being integrated across the entire mush, it can be replaced by the spacial derivative of the diffusion equation (eq. B12 ),
In terms of q h , this equation is also
With eqs ( B26 ) and ( B35 ), the upper boundary condition becomes: ) which we simplify further to
Here,
with β eff an effective compressibility of the crystal-poor rhyolite magma, and
is the volumetric flux due to magma withdrawal. A relationship between m h and P can be derived from eq. ( B28 ). At the top of the mush, assuming pressure continuity across the mush-rhyolite interface (i.e. p = P ) leads to ) and after solving for m h :
The upper boundary condition can be expressed explicitly in terms of m . After combining eq. ( B8 ) with eq. ( B9 ), as well as eq. ( B37 ) with eq. ( B41 ), we find
With eqs ( B12 ), ( B15) and ( B42 ), the entire mush system has been defined. An initial condition is now necessary to solve. For this problem, the mush and rhyolite layer are assumed to be stagnant at the beginning, meaning there are no gradients in pressure or in m . Therefore, the initial condition is
This indicates no flow is occurring in the mush.
The predicted change in porosity allows one to assess whether poroelastic assumptions are holding. The change in porosity, φ, is defined as
where φ f is the final porosity. An expression for φ f can be derived by starting with the base definition of porosity: ) with V f being the final volume of the pores and V being the final total volume of the mush. One can define V as the initial volume plus any change in the mush thickness,
Since volume can be defined by density and mass,
where M 0 is the initial mass of the pore fluid in the mush, which is φ 0 ρ 0 V 0 . Substituting B17 in for M mush and approximating ρ gives
Di viding b y the expression for V gi ves
Since u h is small compared to h , the leading multiplicative term is approximately one. Therefore, the time dependent change in porosity can be approximated as
The mass leaked from the mush has been defined two ways: using the Darcy flux of melt out of the top of the mush (eq. B19 ) and based on the total change in pore fluid content inte grated ov er the mush (eq. B17 ) . Ideally, these values will be equal at every time step, ho wever , numerical error makes this not often the case. The mass balance equation is thus
All previous equations have been in dimensional form. The dimensionality components are
where primes indicate non-dimensional variables. Substituting these into eqs ( B12 ), ( B15 ) and ( B42 ), ). The phenocrysts are dominantly plagioclase, olivine, and clinop yrox ene (with minor Fe-Ti oxide), at relative proportions of 0.64, 0.2 and 0.16, respecti vel y. We can estimate K s using the bulk moduli of each phenocryst phase (Biot 1941 ;Bass 1995 ). These are 72 GPa for An 56 , 130 GPa for olivine and 105 GPa for clinop yrox ene (a value between augite and diopside). We thus estimate a bulk modulus of K s = 87.6 GPa based on an average of these bulk moduli weighted by the phase proportions, which is the limit that assumes all grains deform equally (Hickey & Sabatier 1997 ).
C1.2 Mush porosity, φ 0 , and permeability, k
We define porosity, φ 0 , as the volume fraction of enclave not occupied by phenocrysts and estimate an average value of φ 0 = 0.36 from the volume fraction of crystals within the enclaves (W inslo w et al. 2022 ). From φ 0 and the average maximum crystal dimension, a , we can constrain the mush permeability, k , from empirical relations between porosity and permeability of partially crystallized magmas (Cheadle et al. 2004 ;Hersum et al. 2005 ). Here we follow Hersum et al. ( 2005 ) and use the relation
The range in crystal sizes of the enclaves is approximately 1 mm ≤ a ≤ 2 mm, which results in k = 1.7 × 10 -9 m 2 to 6 . 7 × 10 -9 m 2 .
The shear modulus of the crystal mush, G , is estimated based on the work of Kov á čik ( 2001) as
Here, G 0 is the shear modulus of the material that comprises the solid matrix, φ c is the porosity at which the shear modulus goes to zero, and f is a material dependent exponent. From theory, f is predicted to be 2.1 (Sahimi 1994 ), although in practice it typically falls between 1 and 2 (Kov á čik 2001 ). We therefore use a value of f = 1.5. Fur ther more, we use G 0 = 44 GPa, which is for gabbro (Carmichael 1982 ) and φ c = 0.6, which corresponds to the threshold of crystal network formation (Saar et al. 2001 ;Bachmann & Bergantz 2004 ). Given that φ 0 = 0.36, we obtain G = 11.2 GPa.
The Biot-Willis coefficient, α, describes the relation between the change in the volume of fluid stored within the porous material to the change in the material's volume under constant pore pressure. Generally, α = 1 -K / K s . Note that K will al wa ys be less than K s because pores will reduce the strength of the material. Therefore, for an incompressible solid ( K s → ∞ ) one has α → 1. Alternati vel y, if porosity tends towards zero, α → 0. Here we assume α = 0.7 ± 0.2 which, for reference, falls within the range of sandstone at similar pressures (Song & Renner 2008 ) and porosities (Selvadurai 2021 ).
Although S s has a physical definition based on bulk moduli (see Table 1 and Section C1.6), we will treat it as a dependent parameter.
As such, changes in S s are driven by changes in K r , which is the bulk modulus of the mixture of melt with v olatile bubb les and depends on the bubble fraction within the interstitial melt, φ g , as well as the bubble compressibility (eq. 7 ). Approximating the volatile phase as an ideal gas, its compressibility is given by 1/ P ref . Based on Jay et al. ( 2014 ) and Seropian et al. ( 2021 ), we assume P ref = 140 MPa.
We further assume a bulk modulus for the melt of K melt = 13 GPa (Bass 1995 ). Finally, based on the work of Parmigiani et al. ( 2016 ), the largest bubble fraction that can likely remain trapped in a mush by capillary forces is φ g ≈ 0.15. Using these values in eq. ( 7) gives K r > ∼ 8 . 8 × 10 8 Pa.
The storage capacity of the mush generally defines how much melt is added or released from 'storage' as a consequence of a change in the pressure of the pore fluid. Specifically, the uniaxial specific storage, S s , is thus the volume of melt released, per volume of mush, per unit decrease in pore pressure, assuming a constant vertical stress and negligible horizontal strain. It encapsulates the combined compressibilities of the mush matrix and of the interstitial rhyolite melt with volatile bubbles, and it is a key unknown parameter that we seek to estimate. It is given by
C1.7 Melt viscosity, η, and density, ρ 0
We estimate the viscosity of the Cord ón Caulle rhyolite melt as η ≈ 3300 Pa •s. Our estimate uses the viscosity model of Hui & Zhang ( 2007 ) for the known rhyolite composition (Castro et al. 2013 ), a temperature of 908 • C (Jay et al. 2014 ) and a H 2 O content of 4.6 wt per cent (Liu et al. 2005 ), assuming H 2 O saturation at a pressure of 140 MPa (Jay et al. 2014 ;Seropian et al. 2021 ). We use a melt density of ρ 0 = 2300 kg m -3 (Castro et al. 2013 ).
Based on the work of Delgado et al. ( 2016 ) and Novoa et al. ( 2022 ), most of the Cord ón Caulle graben may be underlain by a magma storage complex (Lara et al. 2006 ). This is also consistent with the surface deformation during and after the eruption and estimates for the area encompassed by magma storage, A , obtained from geodetic inversions range from 20 to 200 km 2 (Delgado et al. 2016 ;Novoa et al. 2022 ). Herein we assume a rectangular sill-like storage reservoir located within the Cord ón Caulle graben and with an areal extent of 60 km 2 (Fig. 1 b) and an aspect ratio of 5 (Delgado et al. 2016 ). The ramifications of assuming different areal footprints of Cord ón Caulle's magmatic system are discussed in Section 6.2.
An important parameter in the sill model (Okada 1992 ), upon which the transfer functions u and w are based, is the depth to the top of the rhy olite lay er, d . We use a value of 6 km based on prior estimates (Delgado et al. 2016 ;Novoa et al. 2022 ).
P. R. Phelps et al .
Neither the mush thickness, h , nor the rhyolite thickness, h r are known. Instead, both are estimated.
We use a shear modulus for the host rock, G c , that ranges from more brittle and fractured rocks (1 GPa) to moderately strong rocks (35 GPa) all of which may be reasonable at Cord ón Caulle (Wendt et al. 2017 ;Delgado et al. 2019 ;Novoa et al. 2022 ). We use these values of G c together with a Poisson's ratio of 0.25 to calculate u and w using the sill model of Okada ( 1992 ) (see Section 4.3.1 for details).
The resultant transfer functions have values of 1.3 × 10 -7 ≤ u ≤ 4.6 × 10 -6 m •Pa -1 and w = 0.3324 for the A = 60 km 2 reservoir.
After an exploratory random sampling of all parameters within a predefined range, we found that for any combination of parameters there is a unique value of both S s and η/ k for which the misfit χ is the smallest. For some combinations of parameters the smallest value of χ found within S s -η/ k space corresponds to the minimum in χ that can be achieved for any combination of parameters. We denote this value as χ min . For other combinations of parameters the smallest value of χ achievable within S s -η/ k space is larger than χ min and will be denoted as χ opt . Figs C1 and C2 illustrate the difference between χ min and χ opt . They are based on grid searches over S s and η/ k in the vicinity of the smallest misfit. The value of χ is indicated by the colour, the lighter the colour the lower the misfit. The centre graphs of both Figs C1 and C2 are identical and serve as the reference case with the smallest value of χ corresponding to χ min . Each of the other graphs in either figure corresponds to a change in one parameter, as labelled on the graph, relative to the reference case. The corresponding smallest value in χ for each case differs and is larger than χ min of the reference case and will be referred to as χ opt .
The values of χ opt , which are obtained from a random sampling of all parameters and minimization over S s and η/ k , can be shown to have well defined functional dependencies. These functional dependencies are expressed in terms of three functions, 1 , 2 and 3 , and are shown in Fig. C3 as graphs of χ opt versus i . The small insets in Fig. C3 correspond to the minimum misfit, χ min , and the optimized misfit, χ opt , cases. Each point shown in Fig. C3 represents one parameter combination that results in a misfit of χ opt , as for example case F in Fig. 5 (b). Some of these have a value of χ opt = χ min , as for example case A in Fig. 5 (b).
It should be noted that the three i encompass the three unknown parameters S s , h and h r . It is therefore possible to leverage these functional relations to find which parameter combinations yield the minimum misfit of χ min = 3.75 cm. Specifically, for a given set of a priori parameters that are either constrained to a specific value or to within a defined range (Table 1 ), we seek to find those combinations of S s , h and h r that have χ min as the smallest misfit in S s -η/ k space. To do so we define three second order polynomials that each fit one of the χ opt versus i trends in Fig. C3 . Taking the v erte x of each polynomial as i = i, min , we then algebraically solve for the three unknowns S s , h and h r where χ opt = χ min = 3.75 cm.
To obtain 1 , we turn to the long-term behaviour of the system. At the (hypothetical) steady state, all gradients in pressure and m will have relaxed to the values P ss and m ss , respectively. To find these values, we start with the pressure boundary condition (eq. B37 ) and inte grate ov er time from 0 to infinity:
The left-hand side will simply become P ss . From mass balance (eq. B51 ), the integral over q h can be replaced with Finally, as the eruptive flux goes to zero at the end of the eruption, the time integral over Q gives the total mass erupted, M ∞ E . Eq. ( C4 ) thus becomes
Taking the relationship between P and m h from eq. ( B41 ), at steady state it becomes m ss = ρ 0 P ss S s (1 -γ ) , (C7) which can be substituted into eq. ( C6) for m ss to get
After simplifying and distributing 1/ , we find
Each term in this equation has units of pressure, thus we define the second term on the right-hand side as the instantaneous eruptive pressure change, P E . To understand why this is a pressure, consider the definition of compressibility:
Recalling that = h r β r , since M ∞ E / ( h r A ) is a change mass per volume it approximates an average change in density of the rhyolite layer. Solving for P ss in eq. ( C9 ) gives
and rearranging for the ratio of P ss / P E shows where 1 comes into the long-term behaviour of the system:
Although the primary purpose of considering the long-term behaviour has been to explain the origin of 1 , we also now have an expression that allows us to predict the long-term pressure change and thus surface deformation.
Downloaded from https://academic.oup.com/gji/article/235/1/610/7209152 by Rice University user on 11 July 2023