Alexandra Morand

Geneviève Brandeis

Stephen Tait

Application of a plate model to reproduce surface deformations observed at Uturuncu volcano, Bolivia

Keywords: Surface deformations, Magma storage, Plate Model, Uturuncu

Introduction

Geological context and motivation

Uturuncu volcano, located in Bolivia, has shown an impressive pattern of deformation whose potential to be an eruption precursor has generated much interest. A circular uplift surrounded by an annular subsidence, both centred on the volcano, were revealed by INSAR data between 1992 and 2017 [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF][START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF]. The 60 km width uplift was first observed by INSAR in 1992 [START_REF] Pritchard | A satellite geodetic survey of large-scale deformation of volcanic centres in the central Andes[END_REF][START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF][START_REF] Gottsmann | 50 years of steady ground deformation in the Altiplano-Puna region of southern Bolivia[END_REF], but, this uplift could go back to 1965 [START_REF] Gottsmann | 50 years of steady ground deformation in the Altiplano-Puna region of southern Bolivia[END_REF]. The subsidence has a much smaller amplitude than the uplift and for that reason was not detected initially, but it now seems probable that they started at the same time [START_REF] Gottsmann | 50 years of steady ground deformation in the Altiplano-Puna region of southern Bolivia[END_REF]. The simultaneity of continuous uplift and subsidence is, at first sight, different to what has been observed at many other systems, and therefore most intriguing. The rate of uplift determined from satellite data was fast between 1992 and 2010 (about 1 cm•yr -1 ), then it has slowed down since 2010 and tended toward 0 since 2017 [START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF]. GPS data acquired between 2010 and 2015 also show a slowing down of surface deformation [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. Uturuncu volcano belongs to the Altiplano Puna Volcanic Complex (APVC), which extends over parts of Bolivia, Argentina and Chile. This volcano and the surrounding area were studied during the PLUTONS project, which started in 2009 and whose synthesis was recently published, making this area well documented [START_REF] Kern | Geochronological imaging of an episodically constructed subvolcanic batholith: U-Pb in zircon chronochemistry of the Altiplano-Puna Volcanic Complex of the Central Andes[END_REF][START_REF] Perkins | Topographic constraints on magma accumulation below the actively up-lifting Uturuncu and Lazufre volcanic centers in the Central Andes[END_REF][START_REF] Shen | A one-dimensional seismic model for Uturuncu volcano, Bolivia, and its impact on full moment tensor inversions[END_REF][START_REF] Kukarina | Focused magmatism beneath Uturuncu volcano, Bolivia: Insights from seismic tomography and deformation modeling[END_REF][START_REF] Farrell | Seismic attenuation, time delays, and raypath bending of teleseisms beneath Uturuncu volcano, Bolivia[END_REF][START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF][START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF][START_REF] Naranjo | Volcanism and tectonism in the southern Central Andes: Tempo, styles, and relationships[END_REF][START_REF] Gottsmann | 50 years of steady ground deformation in the Altiplano-Puna region of southern Bolivia[END_REF][START_REF] Pritchard | Synthesis: PLUTONS: Investigating the relationship between pluton growth and volcanism in the Central Andes[END_REF][START_REF] Mcfarlin | Receiver function analyses of Uturuncu volcano, Bolivia and vicinity[END_REF][START_REF] Lopez | New insights into the magmatic-hydrothermal system and volatile budget 34 of Lastarria volcano, Chile: Integrated results from the[END_REF].

The APVC region has been the location of several cycles of violent eruptions that formed large calderas [START_REF] De Silva | Episodic construction of batholiths: Insights from the spatiotemporal development of an ignimbrite flare-up[END_REF][START_REF] Naranjo | Volcanism and tectonism in the southern Central Andes: Tempo, styles, and relationships[END_REF]. Altogether, 11 calderas, whose diameters range between 10 and 45 km, were formed during the last 25 Ma. They were related to 15 ignimbrite deposits, corresponding to a total magma volume between ≈ 2,950 and ≈ 3,450 km 3 [START_REF] Naranjo | Volcanism and tectonism in the southern Central Andes: Tempo, styles, and relationships[END_REF]. [START_REF] De Silva | Large ignimbrite eruptions and volcano-tectonic depressions in the Central Andes: A thermomechanical perspective[END_REF] estimated that the total volume of erupted magma could be ≥ 15,000 km 3 for the whole APVC. The APVC is linked to a large low velocity zone located at 25 km below sea level [START_REF] Pritchard | Synthesis: PLUTONS: Investigating the relationship between pluton growth and volcanism in the Central Andes[END_REF], whose volume has been estimated at 500,000 km 3 with nearly 25% partial melt [START_REF] Ward | Seismic imaging of the magmatic underpinnings beneath the Altiplano-Puna volcanic complex from the joint inversion of surface wave dispersion and receiver functions[END_REF]. Based on electric resistivity, [START_REF] Laumonier | Giant magmatic water reservoirs at mid-crustal depth inferred from electrical conductivity and the growth of the continental crust[END_REF] showed that this Altiplano Puna Magma Body (APMB) has ultra-hydrated magma, containing about 10 %wt of water. This enormous reservoir could have contributed to the very large wave-length uplift at the scale of the Altiplano [START_REF] Perkins | Surface uplift in the Central Andes driven by growth of the Altiplano Puna Magma Body[END_REF], i.e., much larger than the uplift at Uturuncu which is the focus of the present paper.

One key question raised by this information is whether the spectacular deformation centred on Uturuncu might be precursor to a large silicic eruption similar to those just referred to, involving shallow magma storage. To answer this question requires understanding of the depth and physical mechanism of storage and hence of the stress state in the shallow crust that might lead to instability and eruption. For the ensemble of caldera eruptions referred to above, the average petrologically determined pre-eruptive storage depth of magma is between 4 and 8 km [START_REF] De Silva | Large ignimbrite eruptions and volcano-tectonic depressions in the Central Andes: A thermomechanical perspective[END_REF]. The last eruption of Uturuncu volcano was 250 ka ago. Erupted magmas were mostly dacitic [START_REF] Muir | Experimental constraints on dacite pre-eruptive magma storage conditions beneath uturuncuvolcano[END_REF]. [START_REF] Muir | Experimental constraints on dacite pre-eruptive magma storage conditions beneath uturuncuvolcano[END_REF] showed that the magma-storage depth before eruption was between 3.8 ± 1.9 km below the Altiplano's mean altitude.

Previous models

The presence of the subsidence moat is a key issue because it cannot be explained by a single source in an elastic half space, which has been the standard approach for interpreting volcano deformation for several decades. We also note that the considerable width of the deformation signal, within the latter model framework, implies a deep source. Three more elaborate models have currently been put forward to try to explain the pattern of deformation at Uturuncu volcano. [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF] propose that it is the signal of a rising diapir at the APMB. In their model, the growth of the plume head is responsible for the uplift. The draining of magma from the APMB into the diapir is postulated to cause the subsidence around the uplift. [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] have proposed a system with two sources in an elastic half-space to explain the subsidence around the uplift zone. One source is deep, has negative pressure, and is responsible for the subsidence. The second source is shallower, has positive pressure and is responsible for the uplift. Taking into account the observation of seismic anomalies, [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] considered in their model two sources at 20 km and 70 km depth respectively. To reproduce correctly the observed deformations, the ratio between the volume change at depth and the volume increase of the shallow source is 10:1. This solution requires that the deep source loses a volume ≈10 times greater than that causing pressurization of the shallow source, which is hard to reconcile with mass conservation. We note also that neither of these sources correspond to the much shallower depths of pre-eruptive storage referred to above, albeit that this is not a requirement. Finally, [START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF] suggest the deformation may be due to the transport through the crust of fluid exsolved during mush destabilization in a column-shaped structure extending from APBM. They further suggest that subsequent degassing of these exsolved volatiles leads the deformation to reverse, in order to explain why sediments in the surrounding lakes do not appear to show a permanent tilt [START_REF] Perkins | Topographic constraints on magma accumulation below the actively up-lifting Uturuncu and Lazufre volcanic centers in the Central Andes[END_REF]. We will discuss further these models later, and conclude these preliminary remarks by noting that the previous modelling attempts all require quite a lot of complexity even to capture the first order features of the data and hence that there is room for a fresh approach.

1.3. Motivation for a plate model [START_REF] Kay | Geochemical, isotopic and single crystal 40 Ar/ 39 Ar age constraints on the evolution of the Cerro Galán ignimbrites[END_REF] showed, in the case of Cerro Galán (a similar case further to the south) that the silicic magma was generated by ≈50% of crustal melting and ≈50% from differentiation of more mafic magma. Whether the silicic magma be produced by differentiation of mafic magma or by partial fusion of the crust, the presence of silicic magma implies that the crust is abnormally hot, and hence that the brittle-ductile transition should be shallow. In this region it is estimated to lie between 4.6 and 12.5 km [START_REF] Jay | Shallow seismicity, triggered seismicity, and ambient noise tomography at the long-dormant Uturuncu Volcano, Bolivia[END_REF], based on heat flow measurements and assuming different crustal lithologies [START_REF] Jay | Shallow seismicity, triggered seismicity, and ambient noise tomography at the long-dormant Uturuncu Volcano, Bolivia[END_REF]. Seismic data, however show very little seismicity beneath approximately 5 km [START_REF] Jay | Shallow seismicity, triggered seismicity, and ambient noise tomography at the long-dormant Uturuncu Volcano, Bolivia[END_REF] , and hence this value has been widely accepted in the literature as a plausible depth for the brittleductile transition [START_REF] Hickey | The large-scale surface uplift in the Altiplano-Puna region of Bolivia: A parametric study of source characteristics and crustal rheology using finite element analysis[END_REF][START_REF] Henderson | Decadal volcanic deformation in the central andes volcanic zone revealed by InSAR time series[END_REF][START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. Moreover, it appears to be relatively constant across the Altiplano based on measurements of high heat flow [START_REF] Springer | Heat-flow density across the central Andean subduction zone[END_REF] and associated thermal models of the crust after emplacement of the APMB [START_REF] De Silva | Episodic construction of batholiths: Insights from the spatiotemporal development of an ignimbrite flare-up[END_REF]. Magma storage models based on an elastic half-space are not consistent with these constraints.

In the region of Uturuncu, the crust is thus likely to be composed of an elastic and brittle layer underlain by ductile material. Indeed the layer overlying the BDT can, to first order, be considered to be an elastic plate of essentially infinite lateral extent which is relatively thin compared with the crustal thickness. [START_REF] Pérez-Gussinyé | Effective elastic thickness variations along the Andean margin and their relationship to subduction geometry[END_REF] show the calculation of the elastic thickness along the Andes based on inversion of gravity data. They point out that the elastic thickness is around 20 km below the Altiplano-Puna plateaus. But other authors also show that the implementation of the APMB in the crust, should result in a reduction of the elastic thickness just above it [START_REF] De Silva | Episodic construction of batholiths: Insights from the spatiotemporal development of an ignimbrite flare-up[END_REF][START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF] . This reduction is due to the increase of crust temperature around APMB. The Uturuncu volcano as show by [START_REF] Muir | a. Petrological imaging of an active pluton beneath Cerro Uturuncu, Bolivia[END_REF] is located in the center of the zone above the APMB. The distance between the Uturuncu and the boundary of the zone above APMB are ≈ 60 km to the North, ≈ 160 km to the East, ≈ 100 km to the South and ≈ 130 km to the West. As the total radial deformation in Uturuncu is ≈ 60 km, we admit that the hypothesis of an infinite plate is questionable in the North, but it is good approximation the East, South and West.

The value of ≈5 km for the elastic-plate thickness is not crucial here, indeed inversions detailed below cover a range of 4.5 to 12.5 km. The main point is the geometry of the considered storage, and by extension, that of the deformation source. That this thickness is of the same order as the pre-eruptive storage depth of the last eruption of Uturuncu and numerous large forming-caldera ones in the region, given their ages, may or may not be a coincidence. Nevertheless, these observational constraints may have profound implications for magma storage and motivated us to construct a new model framework which we test in this paper with the Uturuncu data. The deformation pattern (simultaneously involving uplift and subsidence) is reminiscent of that of an elastic plate. We decided to adapt a model from the literature to describe the deformation that would result from the storage of buoyant magma beneath a plate.

Finally, based on gravity data, structures resembling diapirs connected to the APMB were identified by [START_REF] Del Potro | Diapiric ascent of silicic magma beneath the Bolivian Altiplano[END_REF]. [START_REF] Del Potro | Diapiric ascent of silicic magma beneath the Bolivian Altiplano[END_REF] emphasize that dykes would be of a size beneath the resolving power of the method, giving support to a working hypothesis involving diapirs.

Rather than explaining the uplift at Uturuncu as due to growth of a diapir from the APMB, as done by [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF], we explore the idea that the deformation may result from the slow arrival of a diapir of buoyant magma beneath a plate on which it impinges.

Experimental reproduction of subsidence around a central uplift

In this section, we provide a qualitative experimental demonstration that buoyant fluid trapped beneath an elastic plate can produce uplift surrounded by subsidence. We performed lab experiments to reproduce this storage configuration in a simplified way.

Experimental setup

The silicic-magma diapir is expected to rise through the crust at low speed, such that it deforms slowly the ductile part of the crust. Because of this, we consider to first order that the ductile part of the crust does not accumulate a lot of stress due to the diapiric rise, i.e.. we assume that the ductile crust has a fluid rheology. We thus focus on the stage where magma is stored beneath the plate and we ignore the transient phase during which the diapir arrives just under the plate. This transient stage was already treated by [START_REF] Burov | Asymmetric three-dimensional topography over mantle plumes[END_REF].

Experiments were carried out in a Plexiglas square tank of 0.5 m side and 0.5 m high. The setup is illustrated by figure 1. The ductile part of the crust is represented by sugar solution. This is overlain by a 2 cm thick layer of gelatin, which was poured onto the sugar layer and allowed to solidify. This plate of gelatin has elastic rheology and represents the brittle and elastic upper crust. The Brittle Ductile Transition (BDT) is thus represented by the interface between the sugar solution and the gelatin plate. As explained above, we expect that the BDT is somewhat shallower in the zone beneath the volcanic focus, where diapirs arrive carrying excess heat. This thickness variation at the base of the gelatin layer was achieved in the experiments by positioning a Peltier (electric heater) in the centre of the tank, which allowed us to remelt some of the gelatin plate and create a hollow there.

Because the subsidence has a much smaller amplitude than the uplift that it surrounds, the uplift must be large in order for the subsidence to be measurable. To maximize the buoyancy of the stored fluid and thus produce a large uplift, air was injected below the plate with a needle. The injection was performed near the edge of the tank, i.e. distant from the measured zone at the centre, and then the bubbles were moved towards the cavity at the bottom of the elastic plate.

To measure deformation we use a Moiré method, as used by [START_REF] Delorme | Self-similar growth of a bimodal laboratory fan[END_REF] to measure the topography of an alluvial fan. Data acquisition and calculations of deformation were carried out using the software Light3D. Unlike [START_REF] Delorme | Self-similar growth of a bimodal laboratory fan[END_REF], we used the option that measures the topography of the gelatin surface both before and after air injection, and then calculate the resulting deformation by subtraction of the two.

Observed deformation

The deformation of the whole surface, as obtained from Light3D, is shown by figure 2.a. The central uplift (in red) is surrounded by a moat of subsidence (in blue).

The white colour corresponds to null deformation. To obtain a smoothed profile of deformation, many radial profiles were taken between angles Θ 1 and Θ 2 (grey curves on 2.b) and then averaged (red curve on 2.b). The maximum of subsidence is 0.14 ± 0.05 (mm•yr -1 ), figure 2.b. The injection zone is not visible in figure 2.a. The light-red zone near x = 0 m and y = 0.15 m is due to small air bubbles that got stuck and did not reach the hollow in the plate. Small roughness at (x = 0.14 m, y = 0.15 m) and (x = 0.17 m, y = 0.11 m) are due to imperfections of the gelatin surface.

The maximal uplift observed is 1.43 cm whereas the plate thickness is 2 cm. The deformation is thus comparable with the plate thickness, and so this case is beyond the realm of linear elasticity. Nevertheless, this experiment shows that a buoyant fluid beneath an elastic plate can produce deformation qualitatively similar to that observed at Uturuncu volcano. This result encouraged us to adapt a theoretical plate model to describe more quantitatively the deformation at Uturuncu volcano.

Theoretical model

The model described here is based on the one published by [START_REF] Melosh | On the origin of fractures radial to lunar basins[END_REF][START_REF] Melosh | The tectonics of mascon loading[END_REF]. The basic approach is to solve an elastic problem without making a thin plate approximation, which is useful because we do not know a priori whether the latter is applicable. Hence we are not restricted to the case in which the width of deformation must be large compared with the plate thickness.

The mathematical procedure (assuming axial symmetry) underlying this model is described by [START_REF] Ike | Love Stress Function Method for Solving Axisymmetric Elasticity Problems of the Eleastic Half-Space[END_REF]. The general solution of stresses in the plate produced by an external action is determined. Then, boundary conditions allow us to obtain the specific solution. In the most direct method of solution, the mathematical description is transformed from the real to a frequency domain, using a Hankel transform, and the boundary conditions are defined in the frequency domain. As in [START_REF] Melosh | On the origin of fractures radial to lunar basins[END_REF][START_REF] Melosh | The tectonics of mascon loading[END_REF], we assume linear elasticity, and that the plate is isotropic and homogeneous. Unlike in [START_REF] Melosh | On the origin of fractures radial to lunar basins[END_REF][START_REF] Melosh | The tectonics of mascon loading[END_REF], the external action is applied at the base of the plate and, secondly, the elastic slab has a variable thickness. The configuration of storage supposed is shown in figure 3. In the state before the liquid implementation in figure 3.a, there is no negative topography above the shallow-BDT part as we consider here that the plate density and the lower-medium one are equal (ρ c = ρ down ). This consideration does not change the model adaptation given below, which is also valid if ρ c = ρ down . The physical parameters shown in this figure are defined in table 1, which gives all the notations used in this paper.

We start with the general solution given by [START_REF] Melosh | The tectonics of mascon loading[END_REF], for the stresses in the plate. The vertical stress σ zz (r, z), the radial stress σ rr (r, z), the hoop stress σ θ θ (r, z) and shear

ρ up ρ c ; E ; ν ρ down H 0 H r = H 0 -T exp ( r a ) 2 BDT Before liquid implementation a r T a ρ up ρ c ; E ; ν ρ down H 0 H r = H 0 -T exp ( r a ) 2 BDT u z (r, z = 0) After liquid implementation b r ρ m T l a l Crust (elastic) Air (fluid)
Crust (ductile=fluid) Magma (fluid) stresses σ rz (r, z), σ rθ (r, z) and σ zθ (r, z) are given by where r is the radial position on the r axis, z the vertical position along the z axis, k is wave number and ν is Poisson's ratio. J 0 (kr) and J 1 (kr) are Bessel functions of order 0 and order 1 respectively. Finally, A, B, C, D are the arbitrary parameters that will be determined by boundary conditions.

σ zz (r, z) = - ∞ 0 J 0 (kr)[(A + Cz) cosh(kz) + (B + Dz) sinh(kz)]dk , (1)
σ rr (r, z) = ∞ 0 J 0 (kr) - J 1 (kr) kr [(A + Cz) cosh(kz) + (B + Dz) sinh(kz)] + 2 k J 0 (kr) - (1 -ν) kr J 1 (kr) [D cosh(kz) + C sinh(kz)]dk , (2)
σ θ θ (r, z) = ∞ 0 J 1 (kr) kr [(A + Cz) cosh(kz) + (B + Dz) sinh(kz)] + 2 k νJ 0 (kr) + (1 -ν) kr J 1 (kr) [D cosh(kz) + C sinh(kz)]dk , (3)
σ rz (r, z) = ∞ 0 J 1 (kr) A + Cz + D k sinh(kz) + B + Dz + C k cosh(kz) dk , (4)
σ rθ (r, z) = σ zθ (r, z) = 0 , (5)

210

The vertical and radial displacements, respectively u z (r, z) and u r (r, z), are given by

u z (r, z) = - ∞ 0 J 0 (kr) 1 + ν Ek A + Cz + 2ν -1 k D sinh(kz) dk - ∞ 0 J 0 (kr) 1 + ν Ek B + 2ν -1 k C + Dz cosh(kz)dk , (6)
u r (r, z) = - ∞ 0 J 1 (kr) 1 + ν Ek A + Cz + 2(1 -ν) k D cosh(kz)dk - ∞ 0 J 1 (kr) 1 + ν Ek B + 2(1 -ν) k C + Dz sinh(kz)dk , (7)
where E is Young's modulus.

The applied stress σ appl (r) is due to the positive buoyancy of the fluid beneath the plate. We describe this load using a Gaussian function. This is not a fundamental requirement but is used as a convenient mathematical form for the shape of the base of the plate whose Hankel transform is known. Thus,

σ appl (r) = w exp - r 2 a 2 = w a 2 2 ∞ 0 k exp - a 2 k 2 4 J 0 (kr) dk ( 8
)
where a is the half-width of the Gaussian and w the buoyancy of the magma, i.e.

w = (ρ c -ρ m )g ρ c -ρ up ρ down -ρ up T , (9)
where ρ c is crustal density, ρ m the magma density, ρ up the upper medium density, ρ down the lower medium density, g the acceleration due to gravity, and T the maximum source thickness at r = 0.

The shape of the plate bottom and hence of the source is written as

s(r) = ρ c -ρ up ρ down -ρ up T exp - r 2 a 2 . ( 10
)
Buoyancy acts perpendicularly to the lower plate surface s(r). The applied stress σ appl (r) is broken down into two components, a vertical one σ appl zz (r) z and a horizontal one σ appl rz (r) r .

σ appl (r) =   σ appl (r) sin(γ) r σ appl (r) cos(γ) z   , ( 11
)
where γ is the angle between the horizontal and the lower plate surface defined as

γ = tan -1 d s(r) dr . ( 12
)
In the small angle approximation we have sin(γ) = γ, cos(γ) = 1 and tan(γ) = γ.

Equation ( 11) thus becomes

σ appl (r) =   σ appl rz (r) r = σ appl (r) d s(r) dr r σ appl zz (r) z = σ appl (r) z   . ( 13
)
Hence the vertical component of the buoyancy force is

σ appl zz (r) = (ρ c -ρ m )g ρ c -ρ up ρ down -ρ up T exp - r 2 a 2 , σ appl zz (r) = (ρ c -ρ m )g a 2 2 ρ c -ρ up ρ down -ρ up T ∞ 0 k exp - a 2 k 2 4 J 0 (kr) dk , (14)
and the radial one is

σ appl rz (r) = -(ρ c -ρ m )g ρ c -ρ up ρ down -ρ up T 2 2r a 2 exp -2 r 2 a 2 , σ appl rz (r) = - (ρ c -ρ m )g a 2 8 ρ c -ρ up ρ down -ρ up T 2 ∞ 0 k 2 exp - a 2 k 2 8 J 1 (kr) dk . ( 15
)
Using linear elasticity, we calculate displacements and stresses due to each component and, applying the principle of superposition, sum the two solutions to obtain the total displacement, and the total stress.

Vertical components of the applied stress

For the vertical stress component applied at the bottom of the plate ( 14), the corresponding boundary conditions are

σ zz (z = -H r , r, k) = b zz (k) J 0 (kr), (16)
σ rz (z = 0, r, k) = σ rz (z = -H r , r, k) = σ zz (z = 0, r, k) = 0 . ( 17
)
where H r , the variable thickness of the plate, is

H r = H 0 -s(r) , (18)
b zz (k) is a parameter defined to take into account forces that arise in the deformed state and are opposed to that applied to the plate's lower surface z = H r . It is given by

b zz (k) = w a 2 2 k exp -a 2 k 2 4 1 + F(k) . ( 19
)
where

F(k) = 2(1 -ν 2 )(ρ c -ρ up )g Ek sinh(kH r ) cosh(kH r ) + kH r sinh 2 (kH r ) -(kH r ) 2 , (20)
In many plate models, in which a load is applied to the upper plate surface [START_REF] Haxby | Thermal and mechanical evolution of the michigan basin[END_REF][START_REF] Watts | Isostasy and Flexure of the Lithosphere, chapter 3. Isostasy and Flexure of the Lithosphere[END_REF],

the opposing force is the Archimedes force exerted by the lower medium on the flexed plate. In our case, because the applied force bends the plate upward, the opposing force is the weight of the uplifted plate (ρ cρ up )g which appears in F(k).

Injecting ( 16) and ( 17) into equations ( 1) to [START_REF] Del Potro | Diapiric ascent of silicic magma beneath the Bolivian Altiplano[END_REF] gives

A zz = 0 , B zz = b zz (k) sinh(kH r ) + kH r cosh(kH r ) sinh 2 (kH r ) -(kH r ) 2 , C zz = -kB zz , (21)
D zz = -b zz (k) k 2 H r sinh(kH r ) sinh 2 (kH r ) -(kH r ) 2 .

Horizontal components of the applied stress

The radial component applied at the bottom of the plate ( 15) is expressed as a shear stress for which the boundary conditions are

σ rz (z = -H r , r, k) = b rz (k) J 1 (kr), (22)
σ rz (z = 0, r, k) = σ zz (z = -H r , r, k) = σ zz (z = 0, r, k) = 0 , (23)
where b rz (k), like b zz (k), is a parameter that is defined to take into account the force generated in the deformed state that opposes the radial component of stress applied at the base of the plate. The parameter b rz (k) is given by

b rz (k) = - (ρ c -ρ m )g a 2 8 ρ c -ρ up ρ down -ρ up T 2 k 2 exp -a 2 k 2 8 1 + F(k) . (24)
As previously for b zz (k), the 1 + F(k) term in b rz (k) represents the reduction of the applied stress (8) due to the plate's weight which is opposed to the upward push of the fluid.

Injecting ( 22) and ( 23) in equations ( 1) to [START_REF] Del Potro | Diapiric ascent of silicic magma beneath the Bolivian Altiplano[END_REF] gives

A rz = 0 , B rz = -b rz (k) kH r sinh(kH r ) sinh 2 (kH r ) -(kH r ) 2 , C rz = -kB rz , (25)
D rz = -b rz (k)k sinh(kH r ) -kH r cosh(kH r ) sinh 2 (kH r ) -(kH r ) 2 ,
4. Method

Inversion method

Data points for Uturuncu were extracted from figure 4 of [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. This figure shows the mean annual vertical deformation inferred on an inversion of multiple INSAR observations between 1992 and 2010. This corresponds to the maximum uplift rate, because the uplift has slowed down since 2010 [START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF]. These points were used to perform an inversion with the model described above.

Although these data are in fact speed of deformation, we will first treat them as static deformation, and then discuss the non-static aspect of the problem. Data are inverted to find optimal values of a and T , respectively the half-width and the maximum height of the Gaussian-shaped source. The python function scipy. optimize.

curve fit() was used to perform the inversion. This function fits data with non-linear least squares method by adjusting parameters of a given model. Here, our model is the plate deformation model described above and parameters that should be adjusted are the height T of the Gaussian source and its half-width a.

Data half-width is 21.4 km for "negative" r and 23 km for "positive" r. The baseline half-width for data is thus taken as 22.2 ± 0.8 km. During the inversion, parameter a is allowed to range between 1 km and 30 km. The maximum uplift given by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] is 9.1×10 -3 ±0.5×10 -3 m. To reduce the calculation time, the range of possible values for T is also constrained. In well know problem of plate loaded from the top of the plate, the applied stress is compensated with elastic rigidity and the flotability of the bended elastic-part in the lower medium. Depending of the source width, one or the other can be dominant. In this problem, the applied stress come from the bottom, so it has compensated with the elastic rigidity and the weight of the lifted part of the plate.

The maximum of uplift is reached when the elastic rigidity vanishes and the plate was in the case of isostatic equilibrium. The maximum thickness of the source T can be related to the maximum of deformation u z (r = 0, z = 0) by T (ρ cρ m )g = u z (r = 0, z = 0)(ρ cρ up )g. According to value given in table 2 and that the maximum of observed uplift is around 0.01 m, we have T ≈ 0.06 to 0.14. To be sure not to take a too restrictive upper boundary, we take four times the upper value. Inversion bounds for parameter T are taken as 0 m and 0.6 m.

Range of values for input parameters

To be consistent with our model, which sets z = 0 on the plate surface, all depths reported here are relative to the surface of the elastic plate and thus to the mean elevation of the Altiplano. Some authors give depths relative to mean sea level. Where the latter have been used, published depths have been corrected for the mean elevation of the plateau, which is 4000 m in the region of Uturuncu. [START_REF] Muir | Experimental constraints on dacite pre-eruptive magma storage conditions beneath uturuncuvolcano[END_REF].

Based on magnetotelluric and experimental petrology data, [START_REF] Laumonier | Giant magmatic water reservoirs at mid-crustal depth inferred from electrical conductivity and the growth of the continental crust[END_REF] showed that there could be up to 10wt.% dissolved water in the magma present in the APMB. When this water exsolves in the form of bubbles in the magma, the density of the whole magma body stored under the plate may strongly decrease.

[16] estimate the depth of the brittle-ductile transition, H 0 , to be between 4.5 and 12.5 km. They used both the location of the earthquakes at Uturuncu and heat flow data measured on the Altiplano. [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] summarize in their Figure 15 the values of Young's modulus, Poisson's ratio and crustal density that were determined from seismic data.

Densities considered by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] for the upper crust are between 2400 and 2600 kg•m -3 .

However, [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF] uses a density for the crust as high as 2800 kg•m -3 . They also assume a density contrast between the magma and the surrounding crust of 400 kg•m -3 .

Table 2 summarizes the ranges of values between which the physical parameters could vary based on our literature survey. In order to explore all of them, we will use the dimensionless numbers of our problem.

Dimensionless numbers used to test variation of physical parameters

To non-dimensionalize the equations, we use k = a k, dk = adk, r = r/a, dr = dr/a, z = z(r)/H(r), u i (r, z) = u i (r, z)/H(r), σ i j (r, z) = σ i j (r, z)/G, where G is the shear modulus of the plate. The non-dimensionalization of the model equations highlights four dimensionless numbers. In order to represent our results as a function of all physical-parameter values, we focus on two of these. 2 [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF].

3 [START_REF] Jay | Shallow seismicity, triggered seismicity, and ambient noise tomography at the long-dormant Uturuncu Volcano, Bolivia[END_REF][START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF].

4 [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF][START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF].

5 [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF].

The first one is

β 2 = g(ρ c -ρ m )(ρ c -ρ up )T (1 + ν) (ρ down -ρ up )E , ( 26
)
where g is the acceleration due to gravity, ρ c elastic-crust density, ρ m magma density, ρ up the upper medium density, ρ down the lower medium density, T maximum thickness of magma, ν is Poisson's ratio and E is Young's modulus.

320

The second one is

H 0 = H 0 a , (27)
with H 0 the far-field thickness of the plate and a the half-width of the source.

The dimensionless number β 2 represents the ratio of fluid buoyancy to plate rigidity, i.e. large β 2 implies that fluid buoyancy is large compared to plate rigidity, and Based on the physical parameter values given in table 2, we consider a variation of parameter H 0 between 4.5 and 12.5 km and of parameter β 2 T -1 between 3.96×10 -5 and 5.10×10 -4 km -1 . To test the variation of the parameters, several inversions were made by varying the values of E, (ρ cρ m ), and H 0 each time. For each parameter H 0 and β 2 T -1 , 41 points were considered. A total of 41 2 = 1, 681 inversions were thus performed.

Results

Shape of the sources

The source of deformation in our model has a Gaussian shape, described by its maximum height T and half-width a. The values of a and T , which best reproduce the deformations published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] for the 1,681 inversions, are represented respectively on figures 4 and 5.

In these figures, the effect on the inversion of varying the brittle-ductile transition depth, i.e. variation of H 0 , is read horizontally. The effect of varying elastic parameters or magma buoyancy in the crust is read vertically by the variation of the parameter

β 2 T -1 .
Our results show that the thinner the plate and the greater the magma buoyancy (H 0 small and β 2 T -1 large), the more extensive (large a) and thinner (low T) must be the source in order to reproduce the data published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. Conversely, when the plate is thick and rigid, i.e. H 0 large and β 2 T -1 small, the source that reproduces the data is found to be thick and narrow. Grey areas correspond to calculations whose results no longer respect the 0.6 m limit set for the inversion of parameter T . These results Similarly, the thinner the elastic plate is (H 0 decreasing), the larger a will be. The grey zone corresponds to unrealistic results. The black star corresponds to the calculations shown in figure 6.

have therefore been removed from the figures. The black star shown on figures 4 and 5

represents the calculation that led to figure 6. In order to assess, for all 1,681 inversions, the difference between inversion results

Reproduction of the subsidence

Consider first r 0 , the radial position where the surface deformation transitions from uplift to subsidence. Looking at figure 6, it appears that deformation published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] is not symmetrical with respect to r = 0. For negative r, r 0 is equal to 41.0 km, for positive r, it is equal to 48.5 km. Therefore, we set r 0 for the data of 44.7 ± 3.7 km.

On figure 7a, the difference between r 0 given by the numerical inversion and the data, noted dr 0 , is represented as a function of β 2 T -1 and H 0 . The way of reading this figure is the same as described earlier for figures 4 and 5, however, the color scale is different. Values displayed in green tones are an overestimate of r 0 by the code, while purple tones correspond to underestimates of r 0 , as compared to published data. The color scale is centered on white, representing cases where dr 0 = 0, i.e.. where the code reproduces data very well.

By taking into account the error bar on r 0 of ± 3.7 km, figure 7a shows that our code correctly reproduces r 0 for 248 cases. The maximum difference between our calculations and measured deformations is 7.1 km. Generally speaking cases where either the plate is thin (H 0 small) or source buoyancy is high (β 2 T -1 large) best reproduce the data.

Reproduction of maximum uplift

The second criterion tested is the reproduction of maximum uplift u z (r = 0, z = 0), whose value in the data is equal to 91 ± 0.5 (mm•yr -1 ). The difference between this value and the maximum displacement calculated during inversions, noted du max z , is represented as a function of β 2 T -1 and H 0 in figure 7b. The effect of varying the different parameters listed in table 2 is read as in figures 4 and 5. The color scale is the same as in figure 7a.

Figure 7b shows that the maximum difference between code results and observed data is du max z =0.32 (mm•yr -1 ), which is within the 0.5 (mm•yr -1 ) error bar given by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. In the case of maximum uplift, 338 inversions produced realistic results allowing us to reproduce u z (r = 0, z = 0). As previously for r 0 , data are better reproduced for cases where the plate is thin (H 0 small) or where the source buoyancy is high (β 2 T -1 large).

Reproduction of the maximum subsidence

Next we look at the maximum magnitude of subsidence u min z . Subsidence is better constrained for negative r than for positive r (figure 6). So we have selected as a reference the subsidence value in the negative r, taking u min z to be -0.33 ± 0.5 (mm•yr -1 ).

The difference between the value of u min z from the code and that from the data is represented as a function of β 2 T -1 and H 0 in figure 7c. Since u min z is a negative value of displacement, when the code underestimates the subsidence, the difference is positive.

We note that our model tends consistently to underestimate maximum subsidence with respect to the data. Underestimated values are displayed in purple, but the color scale is not centered on 0 as it was in figures 7a and 7b. It is also true that the error bar on the data of maximum subsidence is relatively large. A total number of 338 inversions gave usable results. Minimum deviation du min z , displayed in figure 7c, is 0.19 (mm•yr -1 ). The maximum deviation is 0.26 (mm•yr -1 ), which remains within the 0.5 (mm•yr -1 ) error bar. Unlike r 0 and u z (r = 0, z = 0), the values of u min z are less well reproduced as the plate becomes thinner and the source has more and more buoyancy, i.e. H 0 decreasing and β 2 T -1 increasing. The white area in figure 7c shows where subsidence is best reproduced. The optimum values of a and T that emerge from the inversion are thus to some degree a trade off between reproducing r 0 and u min z .

Ratios

The ratio of maximum subsidence over maximum uplift is noted R1. To assess the difference between our model and observations, we examine the ratio R1 num /R1 data , where R1 num is the ratio obtained from numerical calculations. R1 data , the ratio obtained from data published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF], is -3.7×10 -2 .

Figure 8a shows variation of R1 num /R1 data versus β 2 T -1 and H 0 . The color scale ranges from purple for the highest values to white for the lowest. In figure 8a the closer R1 num /R1 data is to 1, the better R1 is reproduced by the code. Cases where the ratio is best represented correspond to purple areas on the figure. The R1 ratio given by the code lies between 0.21 and 0.41 times the observed ratio, i.e.. it is significantly underestimated by our code. This result underlines that the amplitude of subsidence shown by the data appears to be high with respect to that of uplift.

The ratio of maximum radial displacement over maximum vertical displacement is noted R2. In the case of InSAR data, the displacement is calculated from the satellite line of sight (LOS) [START_REF] Hickey | The large-scale surface uplift in the Altiplano-Puna region of Bolivia: A parametric study of source characteristics and crustal rheology using finite element analysis[END_REF], which is not strictly vertical. Therefore, there is a difference in the position of maximum displacement as viewed on the ascending satellite orbits and that viewed on the descending satellite orbits. R2 can thus be derived from the difference between the two positions of maximum velocity viewed by the satellites [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF].

We assess the difference between our model and the observations by considering the ratio R2 num /R2 data , where R2 num is the ratio resulting from numerical calculations and R2 data is the one from the data published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. The observed R2 data is 0.30 ± 0.15.

Figure 8b shows variation of R2 num /R2 data versus that of different parameters. The color scale is the same as in figure 8a. The ratio determined by the code is between 0.27 and 0.5 times that observed. Thus, R2 ratio is underestimated by our code, signalling that radial displacement is larger at Uturuncu than is predicted by the model. R1 num /R1 data where R1 num is the ratio obtained from numerical calculations, and R1 data the ratio obtained from published data. b, ratio R2 num /R2 data where R2 num is the ratio obtained from numerical calculations, and R2 data the ratio obtained from published data.

Discussion

Reproduction of the shape of the deformation signal

In assessing the data from Uturuncu, and presenting the results from our inversions, we have placed considerable emphasis on the presence of the moat of subsidence surrounding the uplift. This basic feature of the data is not consistent with any model involving a single source within an elastic medium. As referred to above, one suggestion to overcome this problem has been that there may be more than one source. For example, one can fit the data by assuming a small source with positive excess pressure "embedded" in a larger source with a negative pressure differential [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF] . Although there is no obvious interpretation of such a compound source these authors suggest it may be the signal of a diapir starting to rise from the Altiplano Puna Magma Body (APMB). A deep source undergoing deflation in addition to a shallower source undergoing inflation is also in principle viable. However, the considerable breadth of the deformation field at Uturuncu, requires in elastic half-space models that both sources be relatively deep, in fact much deeper than the probable depth of the brittle-ductile transition (BDT). In finite element models, the existence of subisdence surrounding uplift requires a 25 to 65 year flux between two sources in the ductile medium. The pressure changes associated with this flux produce surface deformation. One way out of the conundrum that sources are in a ductile medium is that the rate of pressurisation of a source nominally beneath the BDT could be fast with respect to the relaxation time of the medium. Nevertheless, taking thermal calculations and temperature and relaxation time dependence

given by [START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF], around the APMB and a few kilometers above it, the relaxation time is on the order of a month. This estimate does not support the "fast pressurization" argument, which requires that the flux be large enough to overcome relaxation during 25 years in order to produce deformations observed by INSAR. However, the viscosity of the APMB is uncertain and so relaxation timescales longer than a month cannot be definitively excluded.

The main goal of the present contribution has been to show that a new kind of single source, namely a simple plate with a thickness variation in which buoyant magma is stored, can explain the data well to first order. From the 1,681 inversions, 338 parameter-combinations based on our new plate model produced a viable source , reflecting the non-uniqueness typical of this kind of problem. Although viable sources could be found for a range of plate thicknesses H 0 , for those between 4.5 and 5.5 km, a source can be found for all values of E, ρ c , ρ m , ρ up and ν listed in table 2. For H 0 larger between 5.5 and 10 km it is always possible to reproduce the data for a large combination of parameter values E, ρ c , ρ m , ρ up and ν, as long as the parameter β 2 T -1 lies between 1.4x10 -4 and 5.1x10 -4 . The elastic parameters used for the example of calculation shown in figure 6 are close to those recently published by [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. The limit case of a thick and rigid plate in the lower right corner of figure 4, 5, 7 and 8, which correspond to H 0 large and β 2 T -1 small, does not correctly reproduce the data, as in those cases, the code oversteps the limit value of 0.6 m for T, so we exclude those results.

Taken together, figures 4 to 8 show that our model reproduces well to first order the vertical deformation as long as the plate is relatively thin and/or the buoyancy of the source is strong (i.e. H 0 small and/or β 2 T -1 large). The horizontal distribution of uplift and subsidence as well as the amount of maximum vertical displacement are especially well constrained by the inversions. This is probably the most compelling feature of the type of model we are proposing. The amount of maximum subsidence is less well reproduced by our model than the uplift, but nevertheless remains within the error bars of the observations. This feature accounts for the fact that our model tends to underestimate the value of R1 with respect to the published data. Figure 8b shows that our model also tends to underestimate R2, the ratio between the maxima of radial to vertical displacement. Maximum vertical displacement is well reproduced, indicating that our calculations work less well for radial displacements.

It is worth noting that [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF] have shown that the ratio between the maximum uplift and the maximum radial displacement, R2, is influenced by heterogeneities close to the surface. For example, the presence of a layer of low rigidity diminishes R2. Our model, which involves a homogeneous, isotropic plate, gives a value of R2 less than that observed, but our model is a very simplified one. Although we have not made any detailed assessment of the possibility that a heterogeneity may be present, we suspect that this basic result would carry over into the kind of plate model that we propose.

Thus the inclusion of a sub-layer (within the plate) of different rigidity, or even the fact that the surface (z = 0) is not strictly flat, could potentially account for discrepancies between the observed and calculated values of R1 and R2. Nevertheless, this simplified model achieves the goal of showing that considering a plate deformation could reproduce the simultaneous uplift and subsidence and their radial distribution to within the error bars of the data.

Implications for a potential shallow magma chamber beneath Uturuncu volcano

The "sombrero" pattern of deformation at Uturuncu has led us to put forward a new but relatively simple physical framework for pre-eruptive magma storage that is coherent with key observational constraints: a shallow brittle-ductile transition, as indicated by seismic and heat-flow data, and petrologic evidence of shallow pre-eruptive storage.

While the last eruption at Uturuncu was not recent, the feature of shallow storage prior to large silicic eruptions appears to be a common denominator for the caldera forming eruptions of the APVC, and arguably even more globally for large silicic systems. Our results based on this new paradigm of storage at the base of a "plate" are promising but preliminary and so merit some discussion.

Our model predicts that the magma layer needed to reproduce the deformation at Uturuncu is thin. Depending on various parameter values, a magma layer in the range 0.083 -0.6 m is enough to reproduce the observed annual uplift. In terms of volume this corresponds to between 0.089 km 3 and 0.092 km 3 per year. These flows respectively involve a broad source with a 0.083 m•yr -1 height increase or a narrower one with a 0.6 m•yr -1 height increase. Considering that these deformations have been observed between 1992 and 2017 [START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF], we estimate that a magma layer of maximal thickness from We suggest that the continuous arrival of magma could correspond to the slow arrival of a diapir that is impinging on the elastic upper crust above the BDT. The static calculations of our model represent this process by a series of sources of increasing thickness. [START_REF] Burov | Ascent and emplacement of buoyant magma bodies in brittle-ductile upper crust[END_REF] found in numerical simulations that for a strong density contrast between magma and crust, rates of rise lie in the range 0.015 -0.063 m•yr -1 . For smaller density contrast they calculated rates of rise in the range 0.015 -0.047 m•yr -1 .The rates of rise implied by our model results are in the range 0.083 -0.6 m•yr -1 . The upper bound is clearly larger than the estimates of [START_REF] Burov | Ascent and emplacement of buoyant magma bodies in brittle-ductile upper crust[END_REF] whereas the lower bound is of the same order of magnitude. We obtain the best fit to the data when we assume strong magma buoyancy.

The APMB is thought to be water-rich and therefore its density may decrease strongly at shallow depth due to volatile exsolution. In our inversion results that gave best fits to the data, strong buoyancy was associated with very broad sources and hence to small speeds of uplift, i.e.. close to 0.047 m•yr -1 . Given the uncertainties concerning rheological parameter values involved, the above values represent good agreement.

We conclude that the deformation observed centred on Uturuncu can plausibly be explained by the slow arrival of a diapir beneath an elastic plate overlying a shallow BDT. The diapir feeds a zone whose overall width is approximately 40 km, corresponding to that undergoing uplift. The surrounding subsidence is caused by the response of the elastic plate to the buoyancy force of the magma. We estimate that between 1992 and 2017 a lense of magma approximately ≈2.25 m has accumulated beneath Uturuncu, and between 4.4 m and 31.8 m if the process occurred from 1965 to 2017, although this quantitative result depends on the precise assumption of magma density and hence volatile content. A magma lense this thin is unlikely to be detectable by seismic methods. Detection of a low resistivity zone by magneto-telluric methods or of a negative gravity anomaly may be more promising, but in either case a dense instrument network capable of high resolution would presumably be necessary. This raises the question of whether this amount of magma is representative of the total amount that is stored at shallow depth or whether it is being added to a much larger amount that was previously present. [START_REF] Perkins | Topographic constraints on magma accumulation below the actively up-lifting Uturuncu and Lazufre volcanic centers in the Central Andes[END_REF] showed that sediments in lakes surrounding Uturuncu do not appear to show a record of progressively increasing tilt which would reflect steadily accumulating deformation. [START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF] interpreted this to mean that the deformation is not permanent, and that the deformation may be largely attributed to the presence of gases whose subsequent escape to the surface allows the deformation to reverse. This may be so. Can we reconcile our hypothesis with the absence of permanent tilt in lake sediments? When a diapir rises, it is a buoyancy-driven flow which induces a return flow in the surrounding ductile medium. It has been shown that little surface deformation occurs until the head of the diapir arrives beneath the BDT, this produces the deformation of both the surface and the diapir itself. According to the value of 0.1 to 5 m•yr -1 given by [START_REF] Burov | Ascent and emplacement of buoyant magma bodies in brittle-ductile upper crust[END_REF],

it takes between 1x10 3 to 100x10 3 years to rise through the first 5 or 10 km above APMB which are ductile. The deformation observed at Uturuncu volcano could be the arrival beneath the elastic part of the crust of a few-meter-high new-diapir, wich only produces deformations when it encounters the non-ductile part of the crust. This diapir is not necessarily the first of a new cycle as magma chambers are constructed by multiple inputs of magma. Currently the rate of uplift appears to be tending towards 0 over the last year or two [START_REF] Lau | Variations in the long-term uplift rate due to the Altiplano-Puna magma body observed with Sentinel-1 interferometry[END_REF]. It will be interesting to follow the evolution of deformation in the years to come to establish whether it reverses or is maintained. In the framework of our model, if the deformation is largely due to bubble-free magma, at least a part of the deformation should remain. However if bubbles are making a large contribution to magma buoyancy and these bubbles eventually escape to the surface, some or maybe a large part of the deformation should reverse along the same lines as the model suggested by [START_REF] Gottsmann | Thermomechanical modeling of the Altiplano-Puna deformation anomaly: Multiparameter insights into magma mush reorganization[END_REF].

Our model is a general framework that may potentially be applied more broadly to silicic systems. Although "sombrero" patterns have not been commonly described, this could be partly explained by the great difficulty in identifying and measuring the very weak subsidence. Indeed strong uplift had been observed for several years at Uturuncu before the presence of the subsidence was detected. Nevertheless, there is evidence for "sombrero" patterns in other places. Another example we have identified is the case of the Socorro Magma Body [START_REF] Fialko | Sombrero uplift above the Altiplano-Puna magma body: Evidence of a ballooning mid-crustal diapir[END_REF] in New Mexico (USA). As is the case for Uturuncu, the terrain there is also a desert area which provides unusually good radar coherence over large distances which facilitates measuring the very small amplitude of subsidence around the uplift. In many or most cases less favourable conditions could easily render the subsidence invisible. It would be interesting to study further the Socorro case to see if our model could reproduce the observed deformations. We also hope that our work may stimulate new observational efforts to search for weak subsidence surrounding strong uplifts to see whether the "sombrero" pattern is more widespread than previously realised.

Conclusions

Centred on Uturuncu volcano there is a circular moat of weak subsidence that surrounds a strong uplift. We assess here whether this "sombrero" pattern could be caused by storage of a buoyant fluid at the rheologic contrast of the Brittle-Ductile Transition (BDT), for example due to the slow arrival of a buoyant magma diapir. In short, the deformation bears the signature of an elastic-plate which prevents the buoyant magma erupting at the surface until a large enough volume has accumulated to make the plate unstable.

To test this hypothesis, we have adapted a plate deformation model and use it to invert published data [START_REF] Henderson | Time-dependent deformation of Uturuncu volcano, Bolivia, constrained by GPS and InSAR measurements and implications for source models[END_REF]. To test the whole plausible range of physical parameter values, we performed 1,681 inversions. The results show that our model correctly reproduces three important features of the data : the radial position of the transition between the uplift and the subsidence r 0 , the maximum of uplift u z (r = 0, z = 0) and the maximum of subsidence u min z . However the model does not reproduce as well the maximum of radial displacement and the ratio between u z (r = 0, z = 0) and u min z . This could be due to the simplicity of the model used, which is a simple elastic plate above a fluid medium. Nevertheless, these more subtle features of the shape of the deformation are likely to be reproducible within this basic kind of model framework by including a bit more complexity, such as variation of elastic parameters. Although non-unique, our results most strongly suggest a plate of ≈ 5 km thickness and strong magma buoyancy.

While the amplitude of the subsidence is inevitably small (which may explain why it has not been more commonly reported), this kind of plate model is particularly powerful for explaining datasets for which shallow (< 10 km) storage has occurred but where the deformation signal and/or the caldera formed is very broad. As these features are quite common for large silicic eruptions, we propose that this new kind of model may be more broadly applicable.

• All data come from published data, which are cited in the text.

• No funding information related to this work.

Figure 1 :

 1 Figure 1: Illustration of the experimental setup used to reproduce deformation similar to that observed at Uturuncu. In order to produce a large uplift, air bubbles (i.e. with large buoyancy) were injected with a needle.

Figure 2 :

 2 Figure 2: Experimental demonstration of a subsidence moat around uplift when a buoyant fluid is stored beneath an elastic plate. Injected fluid is here air. a : Deformation of the gelatin surface after air injection. A moat of subsidence (blue) is visible around a central uplift (red). Blue fringes are residual noise. b : Radial profiles taken between angles Θ 1 and Θ 2 (grey curves) and mean deformation profile (red curve). Subsidence located between r = 0.05 m and r = 0.1 m is highlighted in the top right-corner inset.

Figure 3 :

 3 Figure 3: Schematic illustration of the supposed storage. a: Initial shape of the plate before magma arrival. The variation of the BDT depth is due to local temperature variations. Here we consider ρ c = ρ down , otherwise,a small negative topography should be visible above the shallow-BDT part. b: Deformation of the plate after the arrival of magma. Vertical deformation is exaggerated, this produces the visual difference between a and T on the graphic. Physical parameters shown on this figure are defined in table 1.

Notation

  of the Gaussian source (also that of the BDT * shape) r (m) Radial position (r=0 centred on shallowest part of the BDT * ) E (Pa) Young's modulus H 0 (m) Far-field thickness H r (m) Plate thickness at radial distance r T (m) Maximum source thickness (also maximum variation of BDT * depth) β 2 Dimensionless number ρ c (kg•m -3 ) Density of the elastic-part of the crust ρ down (kg•m -3 ) Density of the ductile-part of the crust ρ m (kg•m -3 ) Magma density ρ up (kg•m -3 ) Density of the medium above the plate * BDT : Brittle Ductile Transition.

  325 vice-versa. H 0 quantifies how the load is applied to the bottom of the plate. When H 0 tends towards +∞ the load is applied under the plate at a point source, whereas small values of H 0 represent broadly distributed loads. The parameter β 2 (26) contains T whereas H 0 (27) contains a. In combination parameters a and T define the geometry of the source, and are those we are trying to determine by inversion of surface deformations. It is not possible, therefore, to know a priori the values of β 2 and H 0 in the case of Uturuncu. The results of the inversion will be represented as a function of H 0 = aH 0 and β 2 T -1 .

Figure 4 :

 4 Figure 4: Values of parameter a determined for the 1,681 inversions. The variability of all physical parameters in inputs is represented by the variation of β 2 T -1 and H 0 . The stronger the source buoyancy is with respect to the plate rigidity (β 2 T -1 increasing), the wider is the parameter a that allows to reproduce the data.

Figure 5 :

 5 Figure 5: Values of parameter T determined for the 1,681 inversions.The variability of all physical parameters in inputs is incorporated in the variation of β 2 T -1 and H 0 . The stronger is the source buoyancy withrespect to the plate rigidity (β 2 T -1 increasing), the bigger is the parameter T that allows to reproduce the data. Similarly, the thinner the elastic plate is (H 0 decreasing), the larger T will be. Grey zone corresponds to unrealistic results. Black star corresponds to the calculations shown in figure6.

Figure 6

 6 Figure6shows that the published deformations are very well reproduced in the case of H 0 = 6.1 km and β 2 T -1 = 2.09×10 -4 km -1 . This β 2 T -1 value was obtained with a Young's modulus E of 21.7 GPa, a crustal density ρ c of 2,710 kg•m -3 , a magma buoyancy in the crust (ρ cρ m ) of 355 kg•m -3 , which implies a magma density of 2.355 kg•m -3 . g is equal to 9.81 m•s -1 . The parameter combinations that lead to the curve shown in figure6are not unique in that the same value of β 2 T -1 can be achieved for different combinations of E, ρ c and (ρ cρ m ). The total number of viable sources that have been determined from the 1,681 inversions is 338, i.e. numerous possibilities for parameters E, ρ c , ρ m , ρ up and ν. Generally speaking, sources that reproduce the deformations, tend to involve thin plates and/or large magma buoyancy relative to plate rigidity.

Figure 7 :

 7 Figure 7: Difference of three parameters between published data from [13] and inversion results. Code overestimated values are in green and underestimated ones in blue and purple. a and b : the white zonecorresponds to points where the distribution is perfectly reproduced. c : note that the color scale is not centered on 0. The grey area corresponds to unrealistic results. The black star corresponds to the calculation shown in figure6. a, difference, dr 0 , for the radial position r 0 where transition from uplift to subsidence occurs. b, difference, du max z , in the case of maximum uplift u z (r = 0, z = 0). c, difference, du min z , in the case of maximum subsidence u min z .

Figure 8 :

 8 Figure 8: Comparison of ratios between published data from [13] and inversion results. The grey area corresponds to unrealistic results. The white star corresponds to the calculation shown in figure 6. a, ratio

2 .

 2 08 m to 15.0 m has accumulated under Uturuncu, representing a total volume from 2.2 to 2.3 km 3 . If we take into account that this uplift could have started in 1965 and ended in 2017, the layer maximal-thickness range is 4.32 m to 31.2 m, which would represent a total volume in the range 4.63 to 4.78 km 3 .

Table 1 :

 1 List of notations used in this paper.

Table 2 :

 2 Range of values taken by the physical parameters of the problem in the case of Uturuncu volcano.

	Physical parameters	Notation (units)	Range of values
	Elastic Plate
	Poisson's ratio 1	ν	0.25 to 0.3
	Young's modulus 2	E (Pa)	10×10 9 to 62×10 9
	Maximum thickness 3	H 0 (m)	4.5×10 3 to 12.5×10 3
	Weight 4	ρ c -ρ up (kg•m -3 )	2400 to 2800
	Source of deformation
	Magma buoyancy (in the plate) 5 ρ c -ρ m (kg•m -3 )	200 to 400
	1 [8, 13].

Acknowledgements

We thank Angela Limare, Thierry Rivet and Abdelhak Souilah for their technical support. We had fruitful discussions with Harsha Bhat and Raphaël Grandin.

Comments from Matthew Pritchard and an anonymous reviewer greatly improved the manuscript.

and published data, we used five key criteria:

• Position of the transition from uplift to subsidence: r 0 , which is shown in figure 6.

• Maximum uplift: u z (r = 0, z = 0), which is shown in figure 6.

• Maximum subsidence: u min z , which is shown in figure 6.

375

• Ratio between u z (r = 0, z = 0) and u min z : R1.

• Ratio between maximum radial displacement and maximum vertical displacement: R2.

Competing interests

• No financial conflicts of interests for any author

• No conflict of interest with respect to the results of this paper.


Rajah Nutakki

Richard Röß-Ohlenroth

Dirk Volkmer

Anton Jesche

Hans-Albrecht Krug Von Nidda

Alexander A Tsirlin

Philipp Gegenwart

Lode Pollet

Ludovic D C Jaubert

Frustration on a centred pyrochlore lattice in metal-organic frameworks

Geometric frustration inhibits magnetic systems from ordering, opening a window to unconventional phases of matter. The paradigmatic frustrated lattice in three dimensions to host a spin liquid is the pyrochlore, although there remain few experimental compounds thought to realize such a state. Here we go beyond the pyrochlore via molecular design in the metal-azolate framework [Mn(II)(ta)2], which realizes a closely related centred pyrochlore lattice of Mn-spins with S = 5/2. Despite a Curie-Weiss temperature of -21 K indicating the energy scale of magnetic interactions, [Mn(II)(ta)2] orders at only 430 mK, putting it firmly in the category of highly frustrated magnets. Comparing magnetization and specific heat measurements to numerical results for a minimal Heisenberg model, we predict that this material displays distinct features of a classical spin liquid with a structure factor reflecting Coulomb physics in the presence of charges.

Over the last few decades, the theoretical study of magnetism on geometrically frustrated lattices has proven highly successful in identifying exotic states of matter, ranging from classical spin ice [1,[START_REF]Spin Ice[END_REF] to quantum spin liquids [3][4][5]. In 3 dimensions, the most well-studied frustrated lattice is arguably the pyrochlore lattice of corner-sharing tetrahedra. For the classical nearest neighbor Ising and Heisenberg models on the pyrochlore lattice, the spin liquid ground states have an elegant description in terms of an emergent U (1) gauge field which leads to their characterization as Coulomb spin liquids [6,7], with excitations interacting via an effective Coulomb potential and characteristic pinch point singularities in the spin structure factor. The clearest experimental realization is found in the spin ice compounds Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 [8], where the local orientations of the spins and dipolar interactions introduce additional energetic (rather than purely entropic) Coulomb interactions, which in turn leads to a description of excitations in terms of magnetic monopoles [9]. In the Heisenberg case, recent experiments on the transition metal pyrochlore fluoride, NaCaNi 2 F 7 , find that it is well described by an S = 1 pyrochlore Heisenberg antiferromagnet (PHAF) and thought to realize a Coulomb-like phase [10]. Unfortunately, it is a recurrent theme for Heisenberg materials that further neighbor interactions, anisotropic exchanges or disorder perturbs the spin-liquid physics at low temperatures [START_REF]Spin Ice[END_REF][START_REF] Hallas | Annual Review of Condensed Matter Physics[END_REF][START_REF] Rau | [END_REF][13][14].

Recently, metal-organic frameworks have emerged as a class of materials for realizing strongly mag-netically frustrated systems [15], offering a new avenue to realizing both familiar and novel geometrically frustrated lattices in the lab. Here, our focus is on the metal-azolate frameworks [M(II)(ta) 2 ], where M(II) is a divalent metal ion and H-ta = 1H-1,2,3-triazole, which exhibit a diamond net with vertices made of M-centred (MM 4 ) tetrahedra [Fig. 1], offering the exciting possibility to engineer novel magnetic structures by inserting additional, magnetically active ions into the pyrochlore lattice. However, so far only a handful of works have addressed their magnetic properties [16][17][18][19][20][21]. In this Letter, we explore the potential of metalazolate frameworks in the context of frustrated magnetism by studying [Mn(ta) 2 ], which realizes a centred pyrochlore lattice. Using a combination of specific-heat and magnetic measurements, abinitio calculations, Monte Carlo simulations, exact diagonalization and analytical insights, we predict the existence of a finite temperature regime around T ≈ 2 K where we expect to find the hallmarks of an underlying classical spin liquid in the minimal model. The correlations in this spin liquid can be understood in the Coulomb framework as a dense fluid of charges created by the center spins, reminiscent of the monopole fluid studied in the context of spin ice [22,23]. Minimal model -Adapting the synthesis procedure of Ref. [24], [Mn(ta) 2 ] was prepared as a white powder sample. Rietveld refinements of Xray diffraction data find that [Mn(ta) 2 ] has the cubic symmetry of the F d 3m space group [17] which we confirm through high-resolution X-ray a, Combined ball-and-stick and polyhedra model of cubic [Mn(II)(ta)2] highlighting the positions of divalent octahedrally coordinated Mn ions arranged in a diamond-type lattice (CSD code: HEJQEV). The metal centres differ with respect to their special crystallographic positions and coordination environments: Mn(1) is located on Wyckoff position 8b (site symmetry 43m), coordinated exclusively by the N2 donor atoms of the µ3-bridging triazolate linker; Mn( 2) is found at Wyckoff position 16d (site symmetry 3m), coordinated exclusively by N1 or N3 donor atoms. b, Schematic representation of the centred pyrochlore lattice with magnetic exchange paths for first (J1, Mn(1)-Mn(2), 3.929 Å), second (J2, Mn(2)-Mn(2), 6.416 Å) and third (J3, Mn(1)-Mn(1), 7.858 Å) neighbors. Mn(1) and Mn(2) are respectively labeled centre (orange) and corner (yellow) spins.

powder diffraction at 5 K. Since the 3d valence band of high-spin Mn(II) ions is half filled, we expect its magnetic moment to be isotropic. Previous conductivity measurements together with periodic DFT band structure calculations [25] classified this compound as a wide-bandgap semiconductor with a bandgap ∆ = 3.1 eV ≡ 36 000 K. We therefore consider [Mn(ta) 2 ] an insulator at temperatures T ≪ ∆. To assess the relevance of different exchange pathways [Fig. 1], we performed DFT calculations assuming isotropic exchanges for high-spin 3d 5 Mn(II) ions between first (J 1 ), second (J 2 ) and third (J 3 ) neighbors [see Fig. 1.b and Supplementary Information [26] ]. We obtain that

J DFT 1 ∼ 2 -4 K and γ DFT ≡ J DFT 1 /J DFT 2 ≈ 1.3 -1.65. The fact that J DFT 1 and J DFT 2
are of the same order of magnitude is likely due to the similar exchange pathways, traversing either two or three nitrogen ions respectively along the triazolate ligand. On the other hand, no such pathway is available beyond second neighbors, hence

|J DFT 3 | < 0.01 K ≪ J DFT 1 , J DFT 2
. Thanks to this separation of scales, it suffices to define a minimal model (CPy) on the centred pyrochlore lattice with only first and second neighbor isotropic couplings,

H = J 1 ⟨ij⟩ S i •S j +J 2 ⟨⟨ij⟩⟩ S i •S j -µ 0 H• i S i . (1)
Confirming that [Mn(ta) 2 ] is described by Hamiltonian (1), requires the appropriate theoretical tools. Fortunately, owing to the large magnetic moment S = 5/2 we can resort to classical Monte Carlo simulations, which have proven to be powerful techniques to describe magnetic properties of pyrochlore materials capable of capturing longrange order as well as unconventional correlations of spin liquids [START_REF]Spin Ice[END_REF][27][28][29][30][31]. Continuous-spin models cannot, however, reproduce specific heat at very low temperatures because entropy is ill defined. There, one needs to consider the discrete nature of quantum spins via e.g. exact diagonalization (ED). Since ED for S = 5/2 is particularly costly in computer time and memory, we will restrict calculations to a single tetrahdral unit of five spins, fitted to high-temperature experimental data. Such an approximation shall provide an independent estimate of the energy scales, to be compared with parameters obtained from DFT and Monte Carlo simulations. In the rest of this paper, S i will denote a 3-component classical spin of length |S i | = 5/2, except when analyzing specificheat data where S i is a quantum S = 5/2 spin. Comparison between experiment and theory -To explore the magnetism of [Mn(ta) 2 ], we measured magnetic susceptibility over three decades in temperature.

Fitting the hightemperature regime (T > 200 K) of the susceptibility χ with a Curie-Weiss law [Fig. 2.a], we obtain an effective magnetic moment µ eff = 6.05 µ B per Mn ion and a Curie-Weiss temperature Θ CW = -21 K. These values agree with a previous report [17] (5.8 µ B and -21.9 K respectively) and with the expected size of the magnetic moment g S S(S + 1) = 5.9 µ B which is a stable value for a Mn(II) ion with high-spin 3d 5 electronic configuration. Θ CW indicates sizeable antiferromagnetic interactions. However, our specific heat measurements do not find a transition until a low T c = 0.43 K [Fig. 2.c]. For comparison, the degree of frustration f ≡ |Θ CW | Tc = 49 is of the same order as the one of the celebrated Kitaev materials [13,32] and substantially larger than the one of most rare-earth pyrochlore oxides [START_REF]Spin Ice[END_REF][START_REF] Hallas | Annual Review of Condensed Matter Physics[END_REF][START_REF] Rau | [END_REF]. It means that not only does [Mn(ta) 2 ] offer an unexplored geometry, but one can also expect a sizeable temperature range above T c where frustration is important. Our main result is the agreement between experiments and Monte Carlo simulations for the (i) magnetic susceptibility in Fig. 2.a, and (ii) magnetization curves in Fig. 2.b, using coupling parameters that confirm the DFT estimates, J MC 1 = 2.0 K and γ MC = 1.51 ± 0.15. Further support for this claim is seen in the fit of the specific heat [Fig. 2.c]: The bump at ∼ 4 K is also consistent with finite-temperature exact diagonalization for J ED correlation length beyond a single frustrated unit. The region 0.43 < T ≲ 4 K thus appears appropriate to support a potentially exotic magnetic structure, where frustrated correlations exist, but have not yet been destroyed by long-range order. We study the nature of this regime in the following.

The Centred Pyrochlore Spin Liquid -To show how the CPy model can host a spin liquid, it is convenient to rewrite Hamiltonian (1) as a sum over (centred) tetrahedra, with index α,

H = J 2 2 α |L α | 2 + const, (2)
with

L α = γS α,c + 4 m=1 S α,m , (3)
where c labels the centre spin while the sum over m runs over the four corner spins. The ground state manifold is therefore defined by minimizing L α = |L α | on all units. For γ ≥ 4, this yields a long-range ordered ferrimagnet with 1/3 saturated magnetization where S α,i=1,..,4 = -S α,c . On the other hand, for γ < 4, the ground state is defined by the local constraint

L α = 0, ∀α. (4)
Such an energetic constraint, familiar from the kagome and pyrochlore lattices, is often used to define a classical spin liquid, provided order by disorder does not select an ordered state. A Maxwellian counting argument [33] gives the degrees of freedom D n of these spin liquids. For the kagome, D 3 = 0 (which is famously marginally disordered), whereas D 4 = N u for the pyrochlore [33], where N u is the number of units making up the lattice.

In the CPy model, we find

D 5 = 3 N u for γ ∼ 1.
Similarly, these spin liquids manifest themselves via a number of flat bands F n as the ground state of their excitation spectrum [34]: F 3 = 1 out of 3 bands for kagome, F 4 = 2 out of 4 for pyrochlore, whereas F 5 = 4 flat bands out of 6 for γ < 4 in the CPy model. Therefore, the CPy model hints at a notably strong spin liquid, even more disordered than the extensively studied ones on kagome and pyrochlore. This is supported by our Monte Carlo simulations [35] showing that CPy lies in the middle of a disordered ground state (found for any γ ≲ 3) which indicates that the ordering mechanism in [Mn(ta) 2 ] lies beyond our minimal CPy model. This is often the case in frustrated magnetism [START_REF]Spin Ice[END_REF][START_REF] Hallas | Annual Review of Condensed Matter Physics[END_REF][START_REF] Rau | [END_REF][13], where perturbations ultimately lift the ground-state degeneracy in materials. In this case, the largest perturbation is likely dipolar interactions with a nearest neighbour strength of 270 mK. Adding these to the CPy model, we find ordering at 250mK in MC simulations, where the ordered state is an unsaturated ferrimagnet and corner spins realizing a planar antiferromagnet on each tetrahedron with the remaining spin weight. Further details are provided in the supplementary information. Our simulations show that the addition of dipolar interactions does not significantly modify the properties of the model in the regime 1 < T < 4 K, where we expect the spin liquid to persist.

An emergent charge fluid -To better understand the nature of this spin liquid and its relation to known pyrochlore physics we look at the magnetic correlations, best visualised through the structure factor in reciprocal space, S(q), which we obtained from our MC simulations [Figs. We find that, in the regime 0 < γ < 2, S(q) is characterized by finite width bow ties, which are found to arise from the correlations between spins residing at the vertices of the tetrahedra. The magnetic structure evolves continuously in this regime, with the width of the bow ties increasing with γ. This suggests that these correlations can be understood in the framework of the Coulomb phase on the pyrochlore lattice. Rewriting the constraint, Eq. ( 4), for each of the spin components a ∈ {x, y, z} as

we interpret ρ a α as a pseudoscalar ("magnetic") charge, with a strength parameterized by γ. Hence, the Coulomb description is that of an effective field coupled to charges on the diamond lattice that break the zero-divergence constraint. For small γ, the charge concentration is low, so Debye-Hueckel theory [36,37] can be used to understand the spin correlations. In this regime, there will be entropic screening of the effective field, resulting in a Lorentzian form for its structure factor B a aa (q a , q b = 0, q c = 0) ∝

1 q 2 a + κ 2 , (6)
where a, b, c ∈ {x, y, z} and Debye-Hueckel theory predicts that κ ∝ γ. At some value of γ, the charge concentration becomes large and we will need to account for additional corrections. Remarkably, we find that the prediction from Debye-Hueckel theory holds in MC simulations up to γ ≈ 1. , displays a number of attractive features such as large, isotropic Heisenberg interactions, an unusually large number of magnetically disordered degrees of freedom (in fact, higher than for most known spin liquids), a reasonably large Ramirez frustration ratio, and a structure factor reflecting Coulomb physics in the presence of charges controlled through γ. Experimentally, external pressure may be a viable tool for controlling γ, similar to spin-1/2 frustrated magnets [38,39]. The consequences of this picture for the inelastic spectrum and the nature of excitations are interesting open questions. Inspired by routes taken for the pyrochlores [START_REF] Rau | [END_REF]28,40], ways to build a multitude of exotic phases by taking a selection of degrees of freedom out of the ground state, e.g., via chemical substitution or with magnetic field, can be thought of. Even more broadly, [Mn(ta) 2 ] belongs to a family of metal-azolate frameworks, whose potential for low energy physics remained till recently to a large extent uncharted. These frameworks provide a versatile platform to engineer (quantum) frustrated magnetism on the centred pyrochlore lattice and beyond, such as [Fe(ta) 2 (BF 4 ) x ] with a degree of frustration f ≈ 27 [19], [Cu(ta) 2 ] with Cu(II) dimers at low temperature [18], and [Cr(II/III)(ta) 2 (CF 3 SO 3 ) 0.33 ] with large exchange couplings [21]. While ligand substitution in similar 1,2,3-triazolate-based Mn(II) networks is known to influence their magnetic properties [41], guest molecule loading in the framework pores should further potentiate the adjustment possibilities in such materials. Thus, they offer us a near infinite playground for design and experimental characterization.

FIG. 1 .

 1 FIG.1.a, Combined ball-and-stick and polyhedra model of cubic [Mn(II)(ta)2] highlighting the positions of divalent octahedrally coordinated Mn ions arranged in a diamond-type lattice (CSD code: HEJQEV). The metal centres differ with respect to their special crystallographic positions and coordination environments: Mn(1) is located on Wyckoff position 8b (site symmetry 43m), coordinated exclusively by the N2 donor atoms of the µ3-bridging triazolate linker; Mn(2) is found at Wyckoff position 16d (site symmetry 3m), coordinated exclusively by N1 or N3 donor atoms. b, Schematic representation of the centred pyrochlore lattice with magnetic exchange paths for first (J1, Mn(1)-Mn(2), 3.929 Å), second (J2, Mn(2)-Mn(2), 6.416 Å) and third (J3, Mn(1)-Mn(1), 7.858 Å) neighbors. Mn(1) and Mn(2) are respectively labeled centre (orange) and corner (yellow) spins.

1 = 1 . 4 FIG. 2 . 1 = 2 .

 114212 FIG. 2.Comparison between experiment and theory for [Mn(ta)2]. a. Curie-Weiss fit (blue, for T > 200 K) of the inverse susceptibility χ -1 obtained from magnetization measurements χ = M/H (black) and MC simulations of over 41000 spins (red). Inset, the inverse susceptibilities at low temperatures. The Curie-Weiss fit yields a Landé factor of g = 2.05. b. Magnetization in an external field over a broad temperature range. MC simulations compare well to experimental magnetization measurements at T > 1 K for J MC 1 = 2.0 K, γ MC = 1.51 ± 0.15. c. The specific heat of [Mn(ta)2] displays a broad bump at ∼ 4 K and a sharp peak at 0.43 K. The former is qualitatively reproduced by full diagonalization for J ED 1 = 1.95 K and γ ED ∼ 1.75 up to a rescaling by a factor of 0.8 along the y-axis. In order to show the γ dependence of Hamiltonian (1), simulation curves in a. and c. are supplemented by dotted and dashed lines corresponding to ±10% variations of γ.

5 FIG. 3 .

 53 FIG.3. Spin structure factor from large-scale MC simulations at T = 1.5 K, γ = 1.5 in the [hhl] plane including all spins (a.) and only sites residing at the corners of tetrahedra (b.). Finite width bow ties at q = (±4π, ±4π, 0) are found in both and the structure factor in (b.). is within 1% of that found at T = 0.07 K. c. The width of the bow ties in the structure factor of the effective field parametrized by κ, see Eq. 6, for different values of the coupling ratio γ at T = 0.07 K. A linear fit was performed to the data up to γ = 1.25 (black line).

  = -γS a α,c = ρ a α

Acknowledgments -We are grateful to Nic Shannon for stimulating discussions. D.V. and R.R.-O. are grateful for financial support from DFG grant VO-829/12-2 (DFG Priority Program 1928 "Coordination Networks: Building Blocks for Functional Systems") and SEM micrographs taken by Ralph Freund. A.A.T., A.J., P.G., and H.-A. K.v.N. were supported by the German Research Foundation via TRR80. L.D.C.J. acknowledges financial support from the "Agence Nationale de la Recherche" under Grant No.ANR-18-CE30-0011-01. R.N., L.P. and L.D.C.J. acknowledge financial support from the LMU-Bordeaux Research Cooperation Program. R. N. and L. P. acknowledge support from FP7/ERC Consolidator Grant QSIM-CORR, No. 771891, and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868. Our simulations make use of the ALPSCore library [42].