Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization.

A 150-km-wide ground deformation anomaly in the Altiplano-Puna volcanic complex (APVC) of the Central Andes, with uplift centered on Uturuncu volcano and peripheral subsidence, alludes to complex subsurface stress changes. In particular, the role of a large, geophysically anomalous and partially molten reservoir (the Altiplano-Puna magma body, APMB), located ~20 km beneath the deforming surface, is still poorly understood. To explain the observed spatiotemporal ground deformation pattern, we integrate geophysical and petrological data and develop a numerical model that accounts for a mechanically heterogeneous and viscoelastic crust. Best-fit models imply subsurface stress changes due to the episodic reorganization of an interconnected vertically extended mid-crustal plumbing system composed of the APMB and a domed bulge and column structure. Measured gravity-height gradient data point toward low-density fluid migration as the dominant process behind these stress changes. We calculate a mean annual flux of ~2 × 10 7 m 3 of water-rich andesitic melt and/or magmatic water from the APMB into the bulge and column structure accompanied by modest pressure changes of <0.006 MPa/yr. Two configurations of the column fit the observations equally well: (1) a magmatic (igneous mush) column that extends to a depth of 6 km below sea level and contains trapped volatiles, or (2) a volatile-bearing hybrid column composed of an igneous mush below a solidified and permeable body that extends to sea level. Volatile loss from the bulge and column structure reverses the deformation, and explains the absence of broad (tens of kilometers) and long-term (>100 yr) residual deformation at Uturuncu. Episodic mush reorganization may be a ubiquitous characteristic of the magmatic evolution of the APVC.


INTRODUCTION
The structural and chemical evolution of continental crust is intimately linked to igneous differentiation processes (e.g., Hawkesworth and Kemp, 2006). Both extensive recycling of material by crustal melting and new additions from the mantle contribute to crustal growth (e.g., Patchett et al., 1982).
The Central Andes is a region where crustal thickening has occurred, leading to the formation of ~70-km-thick crust to the east of the Western Cordillera in the Altiplano-Puna region of southern Bolivia and northern Argentina (lat 18°S to 22°S; e.g., Zandt et al., 1994). Satellite remote sensing studies identified a large (>70 km wide) near circular area of surface deformation (Pritchard and Simons, 2002) centered on Cerro Uturuncu (22.270˚S, 67.180˚W), a 6008-m-high Pleistocene dacitic volcano (e.g., Sparks et al., 2008) in the Altiplano-Puna volcanic complex (APVC; de Silva, 1989). The mean maximum deformation rate of ~1 cm/yr has been ongoing from 1992 to 2011 (Fialko and Pearse, 2012;Henderson and Pritchard, 2013).
As part of the central volcanic zone of the Andes, the APVC is interpreted as the volcanic expression of a crustal magmatic system that between 10 and 1 Ma fed episodic and cyclic silicic ignimbrite flare-ups (e.g., Salisbury et al., 2011). The explosive eruption of ~15,000 km 3 of magma across the APVC formed numerous large collapse calderas (de Silva and Gosnold, 2007). Several geophysical methods have detected anomalies under the APVC that are interpreted as the Altiplano-Puna magma body (APMB) (see Table 1 for summary and references).
The deformation anomaly covers a substantial part of the APVC and, although centered on Uturuncu (located ~100 km to the east of the active Holocene volcanic arc), it also includes several other Pleistocene and older volcanic edifices within an area of ~31,500 km 2 (Fig. 1). Several explanations have been proposed for the deformation anomaly, including the pressurization of subsurface cavities (as a proxy for magma reservoirs) located at various crustal depths above, within, and below the APMB, and complex interactions between multiple magmatic systems located at upper to lower crustal levels (Pritchard and Simons, 2002;Sparks et al., 2008;Hickey et al., 2013). Alternatively, Fialko and Pearse (2012) proposed the ascent of a magmatic diapir from the APMB to explain a broad moat of relative ground subsidence surrounding a central area of ground uplift .
Here we integrate and interpret recent results from new geophysical, geodetic, geomorphic, and petrological studies conducted on and around Uturuncu including part of the PLUTONS project (probing Lazufre and Uturuncu together; http:// plutons .science .oregonstate .edu /about; a collaboration core-funded by the U.S. National Science Foundation and the U.K. Natural

Petrology
The piecemeal assembly and eruption of chemically distinct batches of dacite magma from a chemically heterogeneous magma system beneath Uturuncu dominated the eruptive history during its protracted activity from 1050 to 250 ka (Muir et al., 2015(Muir et al., , 2014a. Experiments determined the preeruptive dacite storage pressures at Uturuncu to be 100 ± 50 MPa (Muir et al., 2014b) at ~870 °C (1143 K). Isotopic data from the Uturuncu dacites indicate a mixture of crustal-derived and differentiated mantle-derived melts (Michelfelder et al., 2014). Andesitic and basaltic andesitic enclaves in Uturuncu lavas and in silicic ignimbrites of the APVC (Michelfelder et al., 2014) are interpreted as aliquots of the least evolved resident magmas in the APMB entrained into the volumetrically dominant, silica-rich dacite magmas (Sparks et al., 2008). The mafic magmas have less radiogenic strontium and more radiogenic neodymium isotopes than associated dacites (Michelfelder et al., 2014) and are thought to be Size and geometry of anomalous structures Schilling et al. (2006); Comeau et al. (2015Comeau et al. ( , 2016 Gravimetry A pproximately 12-24 km wide, vertically elongated (8-12 km high), low-density structures rooted at the top of the APMB at ~15 km below sea level Absence of sustained residual long-term deformation Null result for residual ground deformation at time, t > 100 yr Perkins et al. (2016) Note: See text for details on numerical model. APMB-Altiplano-Puna magma body.
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 the heat source for crustal anatexis (Muir et al., 2015). These observations and interpretations are consistent with the hot zone concept of Annen et al. (2006), whereby repeated incursions of hot mantle-derive magmas into the lower or mid-crust drive differentiation by a combination of crystallization and anatexis.
Recent high-pressure electrical conductivity and petrological experiments (Laumonier et al., 2017) provide further insights into magma chemistry in the APMB. They show that the high electrical conductivity (~1 S m -1 ) of the APMB, while maintaining the low melt fractions defined seismically (Ward et al., 2014), can be matched by andesitic melts with high water contents (7-10 wt%). Note that an alternative scenario, wherein the APMB contains dacite melt containing the observed dissolved water content of ≤4 wt%, cannot explain the high electrical conductivity of the APMB as imaged by Comeau et al. (2015). To do so, a substantially higher (by at least a factor of 1.5) melt fraction compared to the seismically constrained fraction would be needed (Ward et al., 2014). Phase equilibrium experiments on andesite enclaves from Uturuncu show that the observed phenocrysts assemblage and characteristic calcic plagioclase composition require pressures in the range 500-700 MPa, magmatic temperatures of 970 °C (1243 K) and dissolved H 2 O contents of 8-10 wt%, consistent with the conductivity studies (Laumonier et al., 2017). Together, the geophysical and petrological data support an APMB composed of 8-30 vol% water-rich andesite melt dispersed in a solid matrix. Mixing of these andesitic melts with crustal melts atop the APMB gave rise to the hybrid dacites at Uturuncu and elsewhere in the APVC (Fig. 4).

Thermal
The extensive magmatism and volcanism over the past 15 m.y. in the Altiplano region is proposed to have left a thermal fingerprint in the upper crust beneath the APVC (de Silva and Gosnold, 2007). The bottom of the present-day seismogenic zone, inferred to be the brittle-ductile transition zone (BDTZ), around Uturuncu volcano is located at depths of 4.5-10 km beneath the ground  surface and is shallowest directly beneath Uturuncu (Sparks et al., 2008;Jay et al., 2012;Keyson and West, 2013). Pleistocene eruptions of Uturuncu produced dacites having depths of last equilibration before eruption (see following) that coincide with the present-day BDTZ (Muir et al., 2014b). An estimate of the current temperature distribution between the APBM and the BTDZ, constrained by petrological, experimental, and geophysical data, is summarized in Figure 4.  show the deformation field around Uturuncu volcano from InSAR (interferometric synthetic aperture radar) observations between 1992 and 2011. A broad moat of ground subsidence is discernible around uplift centered on Uturuncu. Velocities of line of sight (LOS) uplift have been calculated at a maxi-mum of 1.05 ± 0.1 cm/yr between May 1992 and January 2011, whereas LOS subsidence occurred at a maximum of 0.22 cm/yr (Henderson and Pritchard, 2013). The axis-symmetric uplift has a radius of 35-40 km, and peripheral subsidence extends out to 70-75 km. Reoccupation of a nearby leveling line using the global positioning system (GPS) shows that the ground deformation has been active since at least 1965 at rates comparable to the satellite-derived defor ma tion (del Potro et al., 2013b), but continuous GPS data collected near the center of deformation that span 2010-2015 (Appendix Fig. A1 in Appendix 1) indicate that uplift has slowed to a rate of ~0.24 ± 0.19 cm/yr (Blewitt et al., 2016).

Geomorphology
Geomorphological investigations of shorelines and river terraces reveal an absence of sustained residual long-term and long-wavelength deformation at the geodetically observed rates (Perkins et al., 2016), and identify the onset of Uturuncu's present deformation phase as no more than 100 yr ago. These observations indicate that the recent deformation is a transient pulse.

MODELING RATIONALE AND MODEL PARAMETERIZATION
The geophysical and petrological data imply a mechanically heterogeneous crust, including hot and partially molten regions that will modulate subsurface stresses and strains. This inference has major implications for geodetic modeling. Conventional analytical deformation models that assume crustal elasticity and mechanical isotropy are not suited to explain the observed ground deformation over multiple decades. Here we develop a coupled thermal and geomechanical finite-element model that accounts for inelastic crustal behavior using multiphysics capabilities of the commercial finite-element code COMSOL Multiphysics v 5.1 (https:// www .comsol .com /release /5.1).
Our modeling approach is guided by a number of hypotheses and assumptions constrained by the following observations. 1. The observed ground deformation can be explained by dynamics within the geophysically imaged anomalous structures in the middle and upper crust beneath Uturuncu. 2. The APMB mechanically decouples the upper and middle crust (≤30 km depth from surface) from the lower crust. Processes in the lower crust do not contribute directly to the deformation signature of the APVC. 3. Petrological, experimental, and geophysical data constrain the current temperature distribution between the APMB and the BDTZ (Fig. 4). 4. The APMB is composed of water-rich andesite melt dispersed in a solid matrix. 5. Time-dependent inelastic behavior of a hot upper crust modulates subsurface strain. 6. The net deformation over time scales of centuries is zero. Initial and boundary conditions that are informed by the available multiparametric constraints on the subsurface architecture of the APVC are described in the following and illustrated in Figures 5-7. Model parameters are summarized in Table 2.

Solid Mechanics: Initial and Boundary Conditions
We solve for stresses and strains by invoking a 2D axisymmetric model with radial distance 0 km ≤ r ≤ 125 km and depth 0 km ≥ z ≥ -35 km and solving for ~200,000 degrees of freedom using ~24,000 mesh elements with element size between 1.6 and 2670 m (Fig. 5). The top of the modeling domain, at z = 0 m, represents the summit of Uturuncu volcano (at 6008 m above sea level). Thus to convert to height and depth values relative to mean sea level, a value of 6 km needs to be added to the z values. The lateral extent of the modeling domain is constrained by the size of the seismically imaged APMB and the available spatial InSAR coverage and by computational cost. The 250 km width of the axisymmetric model captures the primary features imaged by the static and dynamic geophysical observations while allowing for rapid convergence of the simulation toward a single solution.
The lower boundary of the modeling domain is at z = -35 km (29 km below sea level). We chose this model dimension in line with assumptions 1 and 2. A 2D axisymmetric model is representative of the primary geophysical observations. Although there are some 3D variations, they are of secondary importance.
The upper model boundary is free of displacement constraints, while the lower and right boundaries have roller conditions that allow free lateral, but no vertical, displacement. The left boundary (r = 0) is the symmetry axis. Four different model domains with different mechanical and thermal properties are initially invoked, i.e., volcanic edifice, crust, column, and APMB (Fig. 5). Pressure changes are prescribed to the column and APMB domains by applying arcuate boundary loads, which induce a differential stress field from equilibrium conditions. Mechanical and thermal properties for all domains are given in Table 2.

Volcanic Edifice
The volcanic edifice of Uturuncu is approximated by a trapezoidal shape, with a 120-m-wide top at z = 0 km and a basal width of 9 km at z = -1.5 km. This geometry simulates a prominence of the edifice of 1.5 km above the mean ele-

Crust
We invoke mechanical heterogeneity of the crust by using 3D shear wave velocities (v s ), mineralogical data, and density data to derive elastic constants for subsolidus regions of the model. First we convert the v s data presented in Ward et al. (2014) to p wave velocity data (v p ) following Brocher (2005): v p (km/s) = 0.9409 + 2.0947v s -0.8206v s 2 + 0.2683v s 3 -0.0251v s 4 . (1) Next we calculate Poisson's ratio (n) and density (r) via: and r (kg/m 3 ) = 1661.2v p -472.1v p 2 + 67.1v p 3 -4.3v p 4 + 0.106v p 5 . ( In order to convert the 3D mechanical data for parameterization of the solid mechanics to the 2D axisymmetric models, we use a third-order polynomial fit to the density data calculated in Equation 3, the mineralogically derived densities presented in Lucassen et al. (2001), and the density model reported in Kösters (1999). The resultant fit (Fig. 6A) is r c (kg/m 3 ) = 2425 -0.04 × z -1.9 × 10 -6 z 2 -3.1 × 10 -11 z 3 .
The fit represents a reasonable approximation to the findings by Lucassen et al. (2001) and Kösters (1999) as well as the 1s upper bound of the data from Equation 3.
We obtain the Young's modulus (E) via where v p is in meters/second, and use a third-order fit to the data to parameterize the crustal Young's modulus (E c ) as a function of depth via: E c (Pa) = 3.6 × 10 10 -6.5 × 10 6 × z -322 × z 2 -0.005 × z 3 .
In the absence of relationships between static (on time scales and amplitudes of observed geodetic ground deformation spanning years to millennia) and dynamic (on time scales and amplitudes of seismic deformation) mechanical properties for the Andean crust, we employ the seismically derived dynamic values. The static Young's modulus may be a factor of a few lower than the dynamic modulus (Gudmundsson, 1983). Use of dynamic properties will therefore yield conservative values for pressure transients toward the upper bound of what may be required to explain the data.
APMB. The laterally extensive main body of the APMB has its top at z = -19.5 km (based on magnetotelluric data; Comeau et al., 2016), a radius of 100 km (based on seismic images; Ward et al., 2014), and a thickness of 10-20 km (Ward et al. 2014;Comeau et al., 2016). The density of the APMB at 1243 K and 500 MPa is calculated as 2630 kg/m 3 , accounting for 25 vol% (Ward et al., 2014) of andesitic melt containing 9 wt% of water (density = 2110 kg/m 3 ) within a solid crystal matrix of density 2800 kg/m 3 . The density of hydrous melts at these conditions is calculated using the parameterization of Ochs and Lange (1999); the solid matrix density corresponds to a quartz gabbro calculated after  Table 2. InSAR-interferometric synthetic aperture radar; r-symmetry axis.
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 Hacker and Abers (2004). The mechanical properties of the APMB domain are summarized in Table 2.
Column. The initial shape and size of the vertical column are determined by the gravity and magnetotelluric data (del Potro et al., 2013a;Comeau et al., 2015). As initial conditions (Figs. 5 and 7) we invoke a cylindrical column with a 6 km radius (a value consistent with the crustal resistivity and density contrast images) that is rooted in the APMB and extends to a depth of z = -6 km (i.e., sea level). The top of the column matches the upper bound of reported seismic brittle-ductile transition depth of between 4.5 and 10 km below the Altiplano (Jay et al., 2012;Keyson and West, 2013;see Fig. 6 and following discussion of thermal boundary conditions). The final shape, orientation, and size of the column and APMB are explored in an iterative forward-modeling approach to fit the observed ground deformation data. We attribute a density of 2580 kg/m 3 to the column based on an average density contrast of -120 kg/m 3 (del Potro et al., 2013a) with the crustal density according to Equation 5 between z = -6 and z = -19.5 km. As we show herein, due to thermal constraints and geophysical observations the column is unlikely to contain partial melt except in its deepest part, where the temperature is thought to be above the watersatu rated dacite solidus (~933 K). Comeau et al. (2016) invoked the presence of saline magmatic fluids in the upper parts of the column at depths of approximately sea level to explain its electrical resistivity of 3-7 Ω·m (Fig. 4). Various volumetric combinations of country rock, magma, fluids, and pore space are consistent with the resistivity and bulk density observations (for further discussion see Comeau et al., 2016;del Potro et al., 2013a). Only the deduced bulk density of the column significantly influences model results via the calculation of lithostatic stresses and resultant stress differentials. The combination of different material fractions for given bulk density has little influence on the model outputs.

Body Load
The models account for lithostatic stresses and gravitational loading, which is invoked by a body load applied to all model domains. In order to balance the gravitational forces, we invoke initial stress conditions on all modeling domains using the cyclic stress loading-compensation sequence (Gottsmann and Odbert, 2014). All domains are consequently free of any significant residual stresses and strains prior to applying the boundary conditions. The residual deformations in all domains at the end of this sequence are negligible and ~10 -10 m.

Initial and Boundary Conditions
The following thermal boundary conditions are invoked. The surface is set to 275 K and the lower boundary at z = -35 km has an inward heat flux of 0.09 W/m 2 . Ignoring radiogenic heat production and assuming constant thermal conductivity, a reference geotherm of 30 K/km is thus obtained (Jaupart and Mareschal, 2009). This reference geotherm is maintained in the model only along the right-side boundary, because the temperature distribution in most of the modeling domain it is strongly influenced by the magmatic plumbing system (APMB and column). We invoke a temperature of 1243 K for the APMB, corresponding to the liquidus temperature of andesite inclusions at Uturuncu as determined by petrological observations and experiments (Sparks et al., 2008;Laumonier et al., 2017). This temperature is consistent with the supposition that the APMB comprises wet andesite melt in a solid matrix, overlain by a partially molten layer containing dacite melt (see Fig. 4). While other deeper storage regions from which andesite has been sourced may exist beneath Uturuncu, our assumption is consistent with the available petrological, geochemical, and geophysical constraints (Michelfelder et al., 2014;Muir et al., 2015;Laumonier et al., 2017). The top of the column at z = -6 km marks the upper bound of the BDTZ, for which we set a temperature of T = 387 °C (660 K) (Priestley et al., 2008). This depth corresponds to the upper bound of experimentally determined confining pressures of 100 ± 50 MPa of pre-eruptive  (Muir et al., 2014b). However, the temperature is considerably lower than the pre-eruptive storage temperature of the last erupted dacites of Uturuncu at 870 °C (1143 K; Muir et al., 2014a) due to secular cooling of the shallow magmatic reservoir over the past 250 k.y. (Fig. 4).
To prescribe thermal boundary conditions to the entire subvolcanic plumbing system at Uturuncu, we reconcile the petrological temperature and pressure data shown in Figure 4 by a linear geotherm of 43 K/km between sea level and the top of the APMB. The resultant temperature distribution across all modeling domains is shown in Figure 7A.

Viscoelasticity
The deduced thermal structure beneath the APVC implies that the mechanical response of the middle and upper crust at low strain rates is time dependent. Here we adopt viscoelastic mechanical behavior using a generalized Maxwell model (for details see Del Negro et al., 2009;Hickey and Gottsmann, 2014). This parameterization consists of one elastic branch, where the mechanical response of the medium is governed by its shear modulus, and one viscoelastic branch with a characteristic relaxation time where A is a constant, H is the activation energy, R is the gas constant, and T is the temperature in K (see Table 2). Figure 7B shows the derived viscosity and relaxation time versus temperature relationship for the crustal rheology.

DEFORMATION DATA InSAR Tracks and Representative Profiles of Deformation
Independent analyses of all available stripmap C-band InSAR data for Uturuncu from the ERS-1/2 and Envisat satellites, spanning May 1992 to January 2011, were presented in Fialko and Pearse (2012) and Henderson and Pritchard (2013). However, there are differences in the data analyzed. Both papers analyzed tracks 282, 10 (both beam 2 for Envisat), and 89 (beam 6 from Envisat only), but Fialko and Pearse (2012) included ScanSAR (scanning synthetic aperture radar) mode Envisat data for track 89, giving additional observations, and Henderson and Pritchard (2013) also included track 3 (beam 2 for Envisat). Models here are compared against stacked InSAR data from these two studies using all four tracks: two ascending tracks (89 and 3), and two descend ing tracks (282 and 10). This data set provides maps of average sur-face ground velocities at ~700 × 700 m 2 resolution with estimated uncertainties of ± 0.4 cm/yr (Henderson and Pritchard, 2013). Each velocity measurement is in the radar LOS, which is sensitive predominantly to vertical, but also horizontal, ground motions. The uplift zone and surrounding moat of subsidence extend ~80-100 km from the center of deformation (Fig. 2).
Axisymmetric finite element models produce a single profile of displacement with Uz (vertical) and Ur (radial) components, which can be compared to any profile extracted from InSAR data by projecting the model output into the radar LOS. The projection is determined by satellite heading and incidence angles (Table 3), which change by <1° and 10°, respectively, within a given InSAR track. For a given azimuth, b, from the center of deformation, model output is compared to LOS according to: where the coefficients a, b, and c map vertical, east-west, and north-south displacement components, respectively, into the radar LOS (e.g., Wright et al., 2004). Because of the north-south orbits of InSAR satellites, north-south motions are not well mapped to LOS: for example, <10% of north-south displacement vectors are mapped into the LOS for Envisat beam 2. East-west displace ments are mapped more efficiently (30%-40% for beam 2), but a purely east-west transect does not include the full footprint of surface deformation. We therefore average InSAR measurements both radially and azimuthally from the center of deformation in order to extract representative profiles of the entire west and east sides of deformation. We then treat our model output as being along an east-west profile such that

Vertical and Radial Components from InSAR Inversion
With the four independent viewing geometries from InSAR data, we solve for vertical, east-west, and north-south components of deformation to compare with model output. However, given the insensitivity to north-south displacements, we use profiles from an east-west swath (10 km north and south of the center of deformation) to solve for vertical and radial deformation components. Our derivation of the radial and vertical components is achieved with different combinations of InSAR data, such that results include a range of permissible profiles. Therefore, our model output is compared against these bounds rather than a single specific profile. Fialko and Pearse (2012) and Henderson and Pritchard (2013) interpreted time series of InSAR data as reflecting a linear rate of uplift ~1 cm/yr between 1992 and 2011 (Appendix Fig. A2 in Appendix 1). However, continuous GPS data confirm a much reduced uplift (Appendix Fig. A1  The linear rates determined by InSAR are averages, and may be tem porally aliased because of the infrequent observations. Short-term changes in the defor ma tion rate could have been missed. All of the InSAR tracks have periods without observations spanning 1-2 yr, and the scatter of cumulative deformation about the best-fitting rate may reflect true displacements or simply noise (Appendix Fig. A2 in Appendix 1). The stacked profiles in Figure 2 present additional evidence for a variable rate of deformation. Independently processed ascending data for track 89 (Henderson and Pritchard, 2013;Fialko and Pearse, 2012) span the same time range, and match very well. However, the uplift rate determined from descending data from track 282 (Fig. 3) changes significantly for three separate time ranges. The red profile, extracted from a single interferogram of ERS data from 12 August 1995 to 31 July 2005, has a much higher rate of defor ma tion (1.1 ± 0.1 cm/yr) compared to that in more recent Envisat data in the same viewing geometry (0.7 ± 0. 1 cm/yr, 2003-2011) or the average 1995-2011 rate (0.8 ± 0.1 cm/yr).

Temporal Deformation Constraints from Geodetic Data
In summary, we seek models that produce an average uplift rate of ~1 cm/yr (but that varies depending on the satellite track) over two decades. However, the models are allowed to have significant variation about this trend on annual time scales. With this in mind, plausible models should be capable of producing significant changes in uplift rate over several years, including a significant decline in uplift rate after 2010.

Time-Dependent Pressure Histories
To match the observations that uplift rates in the center of the deformation may vary on annual time scales about an average rate of 1.05 ± 0.1 cm/yr between 1992 and 2011 and 0.24 ± 0.19 cm/yr between 2011 and 2015, we devise a nonlinear source pressurization term over time. This term introduces piecemeal increments of a dimensionless pressure factor a (shown in Fig. 8A) with smoothed onsets using a continuous first derivative function of a with respect to time, t (i.e., ∂/∂t). In combination with an initial boundary pressure load applied to the respective source domains the term then produces a scaled and stepped pressurization history. The pressure factor versus time history (P vt1) A B Figure 8. Invoked pressure factor versus time (P vt ) histories. The pressure factor scales initial pressure loads on the reservoir boundaries as a function of time. (A) P vt1 that represents the history of the best-fit pressure loads in Figure 10 and the best-fit model in Fig 14. (B) Exploratory histories P vt 2 and P vt 3 for reservoir depressurization. shown in Figure 8A includes 4 pressure pulses over 20 yr and reflects the applied history of the best-fit model.
Considering that there is no long-term (>100 yr) record of residual uplift from the Uturuncu deformation pattern (Perkins et al., 2016), we also test whether uplift is reversible in our model; i.e.; whether the ground uplift trigger by pressurization can be reversed by an equal amount of depressurization. In the absence of observational constraints to determine a depressurization history, we simulate depressurization for a period of 30 yr after the end of pressurization (at t = 20 yr) by a either a stepwise depressurization using a pressure factor parameterization (P vt2) or a linear depressurization (P vt3; Fig. 8B).

Model Fitness
To test model fitness we collapse LOS annual ground velocity data from 8 profiles of descending tracks 10 and 282 (Figs. 2 and 3) across the deformation anomaly into the axisymmetric modeling domain using a manual iterative process. Ascending orbital LOS data are not used initially for goodness-of-fit tests because they are noisier, but are compared with the final model predictions, and so serve as an independent test. Similarly, we compare the model predictions to the radial and vertical velocities at the end of the iterative process. We extrapolate the annual descending velocity data linearly over a 20 yr time frame to obtain a data set that covers the time frame of available satellite geodetic data. This way we can test for time dependence in surface strain partitioning by inelastic responses of the crust.

Single Pressure Source-A Cylindrical Vertical Column
The initial numerical model invoked the pressurization of a single rectangular-shaped column (creating a flat-topped cylinder by rotation about its long axis; Figs. 5, 7A, and 9A) by a finite pressure load. This model reproduced the amplitude of observed ground uplift but did not reproduce the observed peripheral ground subsidence. While some subsidence occurred as a result of time-dependent viscous response of the crust to accommodate the stresses induced by the pressurization, the amplitude of this subsidence was only ~2% of the maximum uplift. This result contrasts with the observed LOS subsidence amplitude of 6%-20% of the maximum LOS uplift. Thus pressure transients within a single and rootless column beneath Uturuncu (i.e., not attached or mechanically coupled to a deeper body) cannot explain the observed ground deformation pattern. We hence abandon further exploration of a rootless column.

Amalgamated Bodies
The observed ground subsidence could be induced by a depressurization of the APMB. In that case ground deformation induced by a pressurized column, which is rooted in a flat-topped depressurizing APMB, would be commensurate with the observed pattern, namely an area of central uplift hybrid column m agmatic column Figure 9. Two-dimensional schematic of the iterative shape refinement process for the plumbing system beneath Uturuncu volcano (z is depth beneath Uturuncu summit; r is radial distance; s.a. is symmetry axis). Figure 8 is applied, yielding the bestfit pressure loads shown in Figure 10. Figure 7A.

(D) The model representation of the best-fit hybrid column (outlined by black bold line), with its top at z = -6 km (corresponding to a lithostatic pressure of 150 MPa and a temperature of 660 K). (E) As in D, but for the best-fit magmatic column (outlined by black bold line), with its top at z = -12.3 km (corresponding to a lithostatic pressure of 330 MPa and a water-saturated dacite solidus of 933 K). The temperature color scale is the same as in
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 surrounded by a moat of subsidence. The predicted transition from uplift to subsidence occurs at a distance of ~6 km from the center of uplift and matches approximately the radial dimension of the column. However, it is crucial that the observed moat of subsidence is located at a distance of ~40-45 km from the center of uplift.
Although pressurization of a column with a radius of 40-45 km coupled with a depressurization of the APMB might explain the observations, we have discounted this scenario for lack of any geophysical evidence of a large shallow body with a diameter of ≥80 km immediately beneath Uturuncu.
In keeping with the dimensions of the anomalous potential field vertical structure (6-12 km radius, 8-12 km vertical extent above the mean APMB depth), we then explored the scenario of a conically shaped APMB, with its cusp directly beneath the base of the column, in generating strains that can explain the observed distance between the center of uplift and the moat of subsidence. We found that although a depressurizing conical APMB achieves some subsidence, neither its footprint nor its amplitude matches the observations. However, invoking a bulge-like rise of the upper surface of an otherwise horizontal APMB that connects to the vertical cylinder (Fig. 9B) creates a surface displacement pattern that is broadly consistent with the observations. Although theoretically a number of other geometrical combinations of the three domains (column, bulge, APMB) coupled with distinct pressure histories may explain the deformation pattern, we focused on finding solutions consistent both with geophysical and geodetic observations. Key to matching the transition from ground uplift to subsidence is the horizontal dimension of the bulge. We find that the distance of the onset of ground subsidence from the center of uplift (~40 km) is roughly determined by the semimajor axis of the pressurized bulge (~35 km in best-fit models). The resultant primary shape and dimensions consist of a pressurized vertical column with a radius of 6 km extending from z = -6 to -19.5 km that connects to a bulging upper surface of the APMB with a radius of 35 km and a maximum height of 0.3 km above the mean APMB depth (Table 4). This assembly then connects to a depressurizing APMB (Fig. 9B).

Shape and Mechanical Property Refinement
While the resultant shapes, sizes, and signs of pressure change of the modeled column and APMB are appropriate to induce, in a broad sense, the observed deformation pattern, the near-field uplift pattern out to ~15 km distance is only poorly fitted by a flat-topped column. Following the systematic evaluation of different upper surface morphologies of the column we found that a domed shape with its top at z = -6 km and a radius of 7.5 km at the depth of transition into the bulge (Fig. 9C) gave a much improved fit.
Given the thermal boundary conditions for the column the resultant temperature distribution within the column indicates that most of this domain is at temperatures below the water-saturated dacite solidus (~933 K). Tracking the 933 K isotherm as a proxy for the presence of partial melt within the column reveals that, at the present time, the best-fit models imply that partial melt only occurs below z = -12 km (i.e., 6 km below sea level). To refine the emerging best-fit model we attribute magma reservoir properties to those parts of the dome-shape structure that are above the solidus temperature while the rest of the structure retains the mechanical properties of solid crust.

Pressure Changes and Data Fits
We explore pressure change amplitudes in the domed column, bulge, and APMB to find acceptable fits to the deformation data invoking a stepped pressure history. Although we apply a uniform pressure load with respect to depth on the domed column, bulge, and APMB, the gravitational loading of the model maps this boundary load into depth-dependent pressure changes whereby the pressure change amplitude (excess pressure) is highest at the top of the column and lowest at its base. In addition, we invoke a decrease in the applied pressure (P ) load with radial distance via a quadratic approximation of dP versus r in the same three domains. The lateral pressure gradient is designed such that the pressure change switches sign from positive in the column and bulge to negative where the bulge transitions into a flat-topped Note: APMB-Altiplano-Puna magma body; z-depth; bsl-below sea level. *Calculated from symmetry axis at r = 0 km; multiply value by 2 to obtain full width. † Pvt1-pressure versus time history 1 (see text Fig. 8).
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 surrounding APMB. Effectively, this procedure simulates a vertical and lateral pressure gradient within the plumbing system along which material flow would be possible. As a result the column and bulge undergo net pressurization with respect to lithostatic pressure, while the surrounding APMB undergoes net depressurization.
In addition to the pressurized bulge and depressurized APMB, we find that the deformation data can be reproduced at equal qualities of fit by two end-member models. (1) A pressurized domed column whose vertical extent is defined by the potential field studies with its top extending to the BDTZ at z = -6 km represents a hybrid structure (hereafter termed hybrid column; Fig. 9D) composed of a suprasolidus, crystal-rich magma reservoir in its lower part and an anomalously hot, but subsolidus, body in its central and upper parts. (2) A pressurized domed column whose top is defined by the solidus isotherm at z = -12 km is a magmatic column that largely represents an igneous mush (Fig.  9E) and may denote the uppermost limit of the APMB according to the seismic data of Ward et al. (2014).
The resultant initial pressure change boundary conditions (in Pa) for the entire plumbing system (AMPB, bulge, APMB, and magmatic or hybrid columns) is derived by a residual minimization procedure between modeled and observed ground deformation to give dP = 95,000 -4.3 × r + 4.5 × 10 -5 × r 2 (12) for the case of a magmatic column and dP = 88,800 -4.2 × r + 4.1 × 10 -5 × r 2 (13) for the case of a hybrid column for 0 < r < 70000 m. These initial conditions are then modulated by the time-dependent best-fit pressure factor function (P vt1; Fig. 8) to yield the best-fit time-dependent boundary pressure loads shown in Figure 10. The resultant initial elastic pressure increments then map into surface strains via the time-dependent crustal rheology. Figure 11 displays the overall fit of the LOS deformation prediction from the models to all InSAR data within a 2D axisymmetric representation of LOS defor ma tion as a function of distance from the maximum uplift. Figure 12 shows the resultant model fits along west-east profiles of the four different InSAR tracks. Decomposing the observed annual LOS data into vertical and radial displacement velocities yields a quasi-symmetric pattern of deformation about the area of maximum uplift (Fig. 13). By mirroring the vertical and horizontal components of the predicted deformation field from the best-fit models along an east-west profile we can demonstrate that our models achieve an acceptable match to the observed data within the observational uncertainties (Fig. 13). In the best-fit scenarios the bulge, magmatic column, or hybrid column undergo the following volume increases, respectively: 15 ± 2 × 10 6 m 3 /yr, 6 ± 1 × 10 6 m 3 /yr, and 5 ± 1 × 10 6 m 3 /yr (Table 4). The uncertainties are derived from least-square fits to volume changes derived from model predictions compared to the observed ground deformation data.

Deformation History-Pressurization
Trimonthly snapshots of the simulated maximum ground deformation show a nonlinear evolution over the 20 yr modeling window (Fig. 14). This reflects the response of the ground surface to both initial elastic stressing and subsequent time-dependent stress relaxation. The parameterization of the relaxation time described in Equation 8 coupled with the applied pressure loads shown in Figure 10 provides a time-dependent deformation history that, by linear fitting, is broadly consistent with the observed mean annual surface veloc ity of 1 cm/yr (Appendix Fig. A2 in Appendix 1).

Deformation History-Depressurization
Invoking the depressurization-time histories (P vt2 and P vt3), whereby the reservoir pressure gradually diminishes with time, we find that for both histories the net ground deformation approaches zero within a few years after the end of the pressure perturbation (Fig. 15). That is to show that in our model all residual stresses and strains can be relaxed by the time-dependent mechanical response of the crust over time scales consistent with a transient episode of deformation at Uturuncu (Perkins et al., 2016). Because no observational constraints are available, we do not discuss further the implications of depressurization.

Best-Fit Model and Relation to Existing Models
Our best-fit model requires a direct physical connection between the APMB, the bulge, and the column. In this model, material transport occurs within a hydraulically coupled plumbing system drawing material from the APMB into the bulge and column.
Our best-fit model gives a net volume change of 21 ± 2 × 10 6 m 3 /yr in the bulge and column (Table 4), which is ~35% less than the annual volume change proposed in Sparks et al. (2008). Required mean amplitudes of excess pressurization to explain the available geodetic time series are very small (<0.006 MPa/yr). The pressure changes are significantly lower than those calculated in earlier models (Pritchard and Simons, 2002;Hickey et al., 2013). This difference is due to several factors: in increasing order of significance, (1) the invoked mechanical heterogeneity in physical properties of the crust, (2) the thermally controlled time-dependent stress relaxation, (3) the size and configuration of the reservoirs, and (4) the invoked bulge above the APMB. The pressurized bulge induces a broad area of uplift at pressure amplitudes of ~10% of the excess pressure amplitude within the column.
Compared to the short-term stress variations over the time scales covered by geodetic observations (years to decades) in our model, en masse diapiric ascent of material from the APMB toward the surface as proposed by Fialko and Pearse (2012) provides an unsatisfactory explanation for the time scale of deformation at Uturuncu. The development of a diapir, its ascent, and ultimately its cessation occur on time scales operating over 10 4 -10 5 yr (e.g., Burov et al., 2003) and are inconsistent with the geomorphological evidence (Perkins et al., 2016).

Nature of Subsurface Stressing
Our best-fit models provide an acceptable match to the observed InSAR data from both ascending and descending orbits. The required pressure changes are modest and can be explained by a migration of material at a volu- Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 metric rate of ~16.5 ± 2 × 10 6 m 3 /yr from the APMB into the bulge and column, which expand by 21 ± 2 × 10 6 m 3 /yr.
We envisage that the APMB and bulge are largely composed of igneous mush, here defined as a mixture of incompressible (crystalline framework with interstitial melt) and compressible (volatile) components. We can also envisage layers of melt and magmatic fluid embedded within mush but too small to resolve by tomography. We postulate that the processes that are causing the deformation involve transfer of material from the APMB into the bulge and column.
In principle, three scenarios can induce depressurization of the APMB and pressurization of the bulge and column. In scenario one (S1) incompressible melt and crystals are transferred from the APMB into the bulge and column, while in S2, compressible magmatic fluids are transferred; S3 is a mix of S1 and S2. The temporal mass and/or density changes in the bulge and column where the stresses behind the ground uplift are created differ significantly between these scenarios (Fig. 16A).
In S1 migration of incompressible melt and crystals with a bulk density broadly similar to the bulge and column densities will transfer mass from a deeper to a shallower crustal level and induce a net mass addition to the bulge and column. As a result, their volume will increase while their bulk density re-mains unchanged; long-term net geomorphic uplift would be expected, such as the pre-400 ka doming observed at Lazufre (Perkins et al., 2016). The bulk density of the APMB will increase and its volume will decrease. The volume loss in the APMB will be compensated by the volume increase of the bulge and column.
In S2 migration of magmatic fluids will induce a net density decrease due to the accompanying depressurization upon upward migration and expansion of magmatic fluids within the bulge and column. As in S1, the bulk density of the APMB will increase while undergoing a volume decrease. Due to the increased expansion of compressible gas in the bulge and particularly in the column in response to lower confining pressures compared to the APMB, their combined density decreases and hence volume increase is larger than the density increase and volume decrease in the APMB. In this scenario the overall mass distribution remains broadly unchanged.
In S3 different combinations of mass and density changes may result from the relative proportions of compressible and incompressible components. For the sake of brevity and relevance, we only explore one example, en masse transfer of APMB mush accompanied by volatile exsolution. This process will result in both a mass increase and density decrease in the bulge and column; the volume loss in the APMB will be smaller than the volume addition in the bulge and column. Long-term net uplift is to be expected, similar to S1.
Petrological constraints on the APMB and the mismatch between the derived volume changes in the APMB on the one hand and the bulge and column  Gravity-height gradient data can help discriminate between these two remaining scenarios. Time-lapse microgravimetric data obtained across the deformation anomaly (see Appendix 2 and Appendix Fig. B1) between 2010 and 2013 do not show any significant changes in gravity beyond the free-air gradient effect (Fig. 16B). This indicates that fundamentally, the overall mass distribution below the anomaly remains the same (within the error of the measurements); however, a density decrease must be invoked (Appendix Fig. B2). This broadly matches the conditions proposed in S2 and points toward surface deformation caused predominantly by water-rich magmatic fluid migration in the subsurface (Fig. 17). The modeled volume increase in the (hybrid or magmatic) column and the bulge amounts to ~21 × 10 6 m 3 /yr. With a total volume of 2600 km 3 for the hybrid column and bulge (or 1900 km 3 for the magmatic column and bulge), this annual volume change amounts to ~0.001% of their combined volumes. Most (~75%) of this volume addition occurs in the bulge that translates into the wide footprint of the uplifted area (~70 km in diameter).
However, the largest excess pressure is in the upper part of the column at confining pressures of ~150 MPa for the case of the hybrid column and ~330 MPa for the case of the magmatic column.

Magmatic Column
The requisite volume increase (and density decrease) can be satisfied with a transfer of low-density material (silicate melt and/or H 2 O-rich fluid) from the APMB into the bulge and column (Fig. 17). Two kinds of silicate melt could ema nate from the APMB: dacite of the type erupted predominantly at Uturuncu and thought to form in response to crustal melting and hybridization atop the APMB, or andesite of the variety recorded as inclusions in Uturuncu dacites and proposed by Laumonier et al. (2017) to be the dominant melt composition within the APMB. The key difference between the two kinds of melt is in their dissolved H 2 O content, which in turn determines the depth at which they would exsolve a water-rich magmatic volatile phase. Uturuncu dacites are water undersaturated, containing ≤4 wt% H 2 O (Muir et al., 2014b) and would not exsolve H 2 O until ~100 MPa confining pressure. Conversely, the andesite Figure 14. Simulated evolution of maximum LOS (line of sight) displacement over a period of 20 yr according to pressure history P vt1 (see Fig. 8) for the hybrid column model with a temporal resolution of 3 months. The InSAR LOS (interferometric synthetic aperture radar) data from t10 (see Table 3) are shown for reference including 1σ error bounds. A linear fit to the modeled data (black squares) over the 18 yr time span covered by InSAR observations (circles) yields a mean velocity of 0.9 cm/yr. This modeled mean deformation rate matches the mean observed InSAR LOS deformation rate of 1.05 ± 0.1 cm/yr within the 95% confidence bounds of the linear fit (see also Appendix Fig. A2 in Appendix 1).

Figure 15. Simulated evolution of maximum LOS (line of sight) displacement for the hybrid column model over a period of 100 yr, according to a superposition of pressure histories P vt1
and P vt2 (see Fig. 8). The evolution of LOS displacements for 20 yr ≤ t ≤ 57.5 yr is shown with a temporal resolution of displacements of 7.5 yr. A net zero residual displacement is predicted within a few years of complete depressurization according to pressure history P vt 2. The same result is obtained for the P vt 3, but is not shown for clarity. The predictions for the model of a magmatic column are essentially the same and are not shown.
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 with 8-10 wt% is close to water saturation (Laumonier et al., 2017) and would exsolve H 2 O almost immediately upon leaving the APMB. We assume here that all of the volume loss of the APMB (16.5 ± 2 × 10 6 m 3 /yr) is due to the transfer of andesitic melt with 9 wt% water (Laumonier et al., 2017) from the APMB into the bulge and column. In this context the magmatic column resembles a cupola structure (Cloos, 2001) protruding from a bulging APMB. For an andesite melt density of 2110 kg/m 3 , 3.5 × 10 10 kg/yr of melt is lost from the APMB and migrates into the bulge and column. At conditions relevant for the upper surface of the APMB (550 MPa confining pressure and 1243 K), the density of water is 631 kg/m 3 (Pitzer and Sterner, 1994). At 330 MPa (i.e., the top of the magmatic column) and 933 K the density of water is 672 kg/m 3 , whereas at the same pressure but 1243 K, its density reduces to 571 kg/m 3 (Pitzer and Sterner, 1994). Regardless of these differences in water density the resultant mass addition to the column and bulge induces a volume increase that matches the mass loss from the APMB to within 2%, and is therefore negligible. This is consistent with having predominantly compressible components involved in the material transfer and matches scenario S2. In addition, there is also the possibility that there is downward transfer of crystal-rich mush to compensate for the upward (potentially decoupled) movement of melt and magmatic volatiles beneath the region of overall inflation (Christopher et al., 2015;Sparks and Cashman, 2017). The magmatic column and bulge can thus be interpreted to represent a vertical protrusion of the upper surface of the APMB.
Based on solubility laws (e.g., Gualda and Ghiorso, 2014), Uturuncu andesite melt at 1243 K and 330 MPa is water saturated at 7.4 wt% water, and therefore ~1.6 wt% of the original dissolved water will have exsolved at this pressure, adding a volume of 0.9 × 10 6 m 3 /yr of water to the column. Adding the different volume terms, i.e., volume lost from APMB (16.5 ± 2 × 10 6 m 3 /yr) + volume of 1.6 wt% exsolved water (0.9 × 10 6 m 3 /yr) gives a mini mum volume addition of 17.4 ± 2 × 10 6 m 3 /yr to the column and bulge. This increase is between 73% and 92% of the modeled average volume increase (21 ± 2 × 10 6 m 3 /yr) in the bulge and column and at the lower bound of error estimates (±10%) of these annual volume changes. Alternatively, if all water is exsolved as the melt solidifies completely a total of 9 wt% of water is released; this adds a volume of 4.8 × 10 6 m 3 /yr to the column and bulge. The total volume change then becomes 21.3 ± 2 × 10 6 m 3 /yr and matches the predicted volume increase (21 ± 2 × 10 6 m 3 /yr) in the column and bulge (Table 4).
These calculations provide very good approximations to the modeled source volume changes and fully explain the volume ratio of 1.3:1 between bulge and column structure inflation and APMB deflation. Due to the different source depths of depressurization (deeper) and pressurization (shallower), this ratio maps into the much larger ratio between the surface volume changes: the volume increase in the area of uplift is a factor of 3-4 larger than the volume of surface subsidence (Fig. 2). Given a compressibility of ~5 × 10 -11 Pa -1 of hydrous andesite melt with 5% H 2 O (i.e., median value) in the absence of a free gas phase (Ochs and Lange, 1997), the volume change by decompression by 220 MPa (between the top of the APMB and the top of the magmatic column) is between 2 and 3 × 10 5 m 3 /yr and can be ignored for these primary calculations.

Hybrid Column
In a hybrid column (composed of a partially molten lower part and a low-density solidified body in its central and upper parts) pressure increase could be caused by the intrusion of water-saturated andesite melt from the APMB that releases all of its water by the time it reaches 330 MPa confining pressure and 930 K, and migration of the liberated water-rich fluids toward shallower crustal  Huppert and Sparks, 2016) indicate that fluid pressures exceed rock strengths at depths of ~10%-30% of their depth of exsolution. Translated to our model, one would therefore expect to see seismicity caused by fluid expansion in the depth range of 4-6 km beneath the surface (i.e., at approximately sea level). This is broadly consistent with the observations at Uturuncu, where seismicity is narrowly concentrated at depths between the surface and sea level immediately beneath the volcano (Jay et al., 2012). This fluid transfer can also explain the exso lu tion and accumulation of magmatic brines at P < 150 MPa (see Driesner and Heinrich, 2007;Weis, 2015) at the top of the hybrid column (Fig. 4).

Fluid Transfer Regimes
The transfer of fluids from large magmatic systems to the shallow crust and Earth's surface involves a very wide range of time scales. At one end of the spectrum, continuous volatile release can occur for hundreds of thousands of years to millions of years during the emplacement of plutons and associated development of magma chambers. Degassing can continue in the fumarolic stage of long dormant volcanoes that remain connected to still active magmatic systems. Greatly enhanced rates of degassing though can occur in much shorter periods (typically months to decades) in association with eruptions and with periods of volcanic unrest induced by magma system processes (Sheldrake et al., 2016). Uturuncu can thus be characterized as a long dormant volcano in its fumarolic stage but undergoing enhanced degassing, as evidenced by the deformation. Our modeling demonstrates the importance of magmatic fluid transport in explaining the observed deformation, shallow seismicity, and associated geophysi cal features such as accumulation of brine at the top of the column.
Transport of magmatic fluids beneath Uturuncu can be divided into two permeable flow regimes; i.e., within the magmatic mush as magmatic fluids are transferred from the AMPB to the bulge and column structure, and within the subsolidus column of rock. Recent observations of degassing and deformation at volcanoes demonstrate that transport can be rapid in igneous mush systems. For example, episodes of magmatic fluid transport from depths of at least 10 km are illustrated by volumetric strain data and SO 2 gas fluxes associated with eruptive activity at the andesitic Soufrière Hills volcano on Montserrat (Gottsmann et al., 2011;Hautmann et al., 2014;Christopher et al., 2015). Transport speeds of several meters per second are indicated within the Soufrière Hills volcano mush region and are consistent with transport through transient fracture systems. Such a rapid fluid transport regime may also be envisaged for the suprasolidus parts of the bulge and column structure be- Figure 17. Schematic synthesis of the model to explain the observed central uplift and peripheral subsidence at Uturuncu volcano. Water-rich magmatic fluids ± andesitic melt are transferred from the Altiplano-Puna magma body (APMB) into a three-phase bulge and column structure containing water-rich fluids, plutonic (solid) residue, and partial melt. Equally good fits to the deformation data are obtained by a depressurization of the APMB that induces peripheral ground subsidence and pressurization of the bulge and column structure that induces broad uplift. Pressurization of either (1) a smaller magmatic (igneous mush) column (outlined by the broken line) formed by melt migration and trapping of exsolved water-rich fluids at above solidus temperature and pressure or (2) a larger hybrid column (marked by solid line) formed by an igneous mush below a subsolidus, altered and fractured plutonic complex through which exsolved fluids migrate toward the surface matches the uplift pattern. Fluids released from the hybrid column may contribute to brine accumulation (Comeau et al., 2015; indicated by wavy area) and seismicity (represented by white stars) above the brittle-ductile transition zone located at approximately sea level (bsl-below sea level) immediately below Uturuncu (Jay et al., 2012). Flow paths of the different components are indicated for illustration by arrows (solid arrows: fluids ± melt; broken arrows: fluids only). Note that the thickness of the bulge is accentuated for clarity and is not to scale. Its maximum thickness in the model is 0.3 km. Anatexis of crustal rocks (shown in brown) above the APMB is indicated by lighter shading.
Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 neath Uturuncu over time scales of tens of minutes to hours. In contrast, in the subsolidus parts of the hybrid column we envisage transport of fluids through a permeable network of interconnected pores and fractures. Damage zones related to the repeated emplacement of dikes between magma chambers and volcanoes and the transient passage of hydrothermal fluids typically enhance permeability by one or two orders of magnitude from typical crustal values in brittle environments (Manning and Ingebritsen, 1999;Ingebritsen et al., 2010;Weis, 2015). In hot but still subsolidus ductile rocks, dynamic permeability develops under conditions of overpressurized fluid flow (Weis, 2015). For typical permeabilities of between 10 -11 and 10 -14 m 2 , transport speeds are between 10 -5 and 10 -2 m/s (Okumura et al., 2009), giving time scales of decades to a few days for fluid ascent through the 6-km-thick subsolidus part of the hybrid column.

IMPLICATIONS AND CONCLUSIONS
The transfer of magmatic fluids and melt from the APMB toward a column-like structure at shallower crustal levels provides a plausible mechanism to explain the evolution of ground deformation at Uturuncu. In our favored model, melt and water-rich fluid are sourced from wet andesitic melts within the APMB at water contents as high as 9 wt% (Laumonier et al., 2017) and transported into a three-phase bulge and column structure containing water-rich fluids, plutonic (solid) residue, and partial melt. The involvement of magmatic fluids explains the mismatch between the volumes of surface uplift and surface subsidence. Subsurface stresses and subsequent surface strains are created either within a magmatic (igneous mush) bulge and column structure by melt migration and trapping of exsolved water-rich fluids, or within a hybrid bulge and column structure by the migration of exsolved fluids from the igneous much into an overlying permeable and altered plutonic complex. The decadal time scale of surface deformation, lack of long-term residual surface deformation, and melt and fluid transfer from mid-crustal depths within the AMPB into the upper crust favor a rapid process of igneous mush organization (Sparks and Cashman, 2017) as the cause for the large deformation anomaly (Fig. 17).
Our findings have the following implications.
1. The upper crustal anomalous body imaged by the potential field studies of del Potro et al. (2013a) and Comeau et al. (2015) beneath Uturuncu could represent a hybrid column with water-rich fluid-bearing fractured and altered plutonic rocks above an igneous mush protruding from the APMB. 2. The presence of several of these geophysically anomalous bodies beneath the APVC (del Potro et al., 2013a;Comeau et al., 2015) suggests that these structures are fundamental for channeling material from the APMB to shallower crustal levels, and may fuel volcanism. The destabilization of segregated volatile-rich melts within the APMB and their subsequent upward migration and fluid exsolution through the anomalous vertical bodies may be a ubiquitous feature of the magmatic and volcanic evolution of the APVC. If this were true one would expect to see a complex temporal evolution of subsurface stresses and surface strains in the APVC over time scales of hundreds or thousands of years. In this scenario, the APVC may undergo complex spatiotemporal surface deformation by the periodic pressurization and depressurization of these vertical bodies. 3. The temporal evolution of the modeled deformation appears to be nonlinear over time scales of months as a result of stepped pressure increments and time-dependent stress relaxation. However, available episodic satellite or ground-based geodetic data lack the temporal resolu tion to resolve such behavior. If periodic material transfer through the crust is responsible for the observed ground deformation, it is imperative to obtain spatially extensive routine geodetic time series at observation frequencies much higher than hitherto available, at least on a monthly basis. With the launch of the European radar imaging satellite Sentinel 1B in 2016, these capabilities might be achievable. 4. Fluid loss from the hybrid column and associated phase separation into vapor and liquid within a shallow-seated, brine-bearing hydrothermal reservoir and concurrent depressurization of the column will reverse the broad deformation pattern and can explain the absence of long-term (>100 yr) and long-wavelength residual deformation around Uturuncu. Periods of ground subsidence and presumed depressurization have been observed at several volcanoes, including Cerro Blanco in the Central Andes (e.g., Henderson and Pritchard, 2013), and could provide clues to what might happen at Uturuncu in the future. 5. Three-dimensional variations in the crustal structure as well as depthdependent regional stresses from the convergence of the Pacific and South American plates and left-lateral transtension above the APMB (Riller et al., 2001) may have an influence on subsurface stress generation and resultant spatiotemporal strain partitioning (see Hickey et al., 2016), which we cannot capture in our axisymmetric model. The first steps to obtain insights into the 3D structure of the middle and upper crust below the APVC and APMB have been made and need expanding. These models in combination with surface deformation data at higher spatial and temporal coverage than currently available will inform future 3D thermomechanical models of the crustal dynamics of the Altiplano-Puna.

APPENDIX 1. InSAR AND GPS TIME SERIES APPENDIX 2. TIME-LAPSE MICRO-GRAVITY SURVEYS
A microgravity network was installed around Uturuncu volcano in March 2010 that consisted of 24 benchmarks that were surveyed relative to a base station (UBAS) located to the west of Uturuncu near the Laguna Colorada. Standard surveying techniques (see Battaglia et al., 2008) were applied, including the establishment of control points as part of individual survey loops to better correct for instrument drift and possible tares.
A Scintrex CG5 spring gravimeter was used for all surveys. Calibration of the instrument was performed before each survey to ensure comparability of results. GPS (global positioning system) receivers were copositioned with each gravimeter benchmark to obtain precise location data during each survey. A total of 5 benchmarks (BLND, KISS, LDZP, MTLC, and OFSP) were occupied during 4 (1 initial and 3  Gottsmann et al. | Mush reorganization at Uturuncu GEOSPHERE | Volume 13 | Number 4 (DEMO, FOOF) were surveyed 5 times (1 initial and 4 repeats) during the same period. A further 17 benchmarks were either installed after March 2010 or reoccupied only between once and twice between 2010 and 2013 and are not considered here. Individual gravity data precision by repeat readings of the same benchmark during individual survey loops was within ±5 and ±15 µGal. Higher uncertainty was usually obtained at benchmarks that were difficult to access due to rough terrain.
Here we only report data from those benchmarks covering the longest survey period (Appendix Fig. B1) along a 60-km-long west-east transect. We use the InSAR LOS (interferometric synthetic aperture radar line of sight) data (see text) to correct for the associated ground deformation. Base station UBAS is located in the moat of subsidence around Uturuncu, where an average LOS subsidence rate of 0.3 cm/yr is detected. We account for this subsidence in the calculation of mean annual elevation change between the base and all other benchmarks.
Calculating the free-air gravity change for the annual minimum and maximum LOS displacements (~90% of vertical displacement) relative to UBAS of 0.3 and 1.2 cm/yr, respectively, yields -4 and -15 mGal/yr. We use the theoretical free-air gravity gradient of -308.6 mGal/m to derive the two gradients. Both gradients are shown in Figure B2 as proxies for gravity changes induced solely by ground deformation, the implication being that observed gravity changes following those trends can be attributed to a change in the subsurface density distribution without a significant mass change.
Most temporal gravity variations observed during the survey period broadly follow the predicted gravity change along the free-air gravity gradient. This indicates that the uplift is predominantly caused by a decrease in the subsurface density while the subsurface mass distribution remains unchanged.