Thomas Reverso
email: thomas.reverso@univ-savoie.fr

Jean Vandemeulebrouck

François Jouanne

Virginie Pinel

Thierry Villemin

Erik Sturkell

Pascale Bascou

A two-magma chamber model as a source of deformation at Grímsvötn Volcano, Iceland

Key Points

Introduction

For basaltic volcanoes, intereruptive ground deformation is often characterized by inflation due to replenishment of a shallow reservoir with magma that was drawn up during the last eruption or during the last diking event, such as for the inflation-deflation cycles during the Krafla rifting episode [START_REF] Björnsson | Rifting of a plate boundary in North Iceland 1975-1978[END_REF][START_REF] Tryggvason | Optical Levelling Tilt Stations in the Vicinity of Krafla and the Krafla Fissure Swarm: Observations[END_REF]. Usually, the inflation rate is greater just after eruptive events and then it decreases with time [START_REF] Lu | Magma supply dynamics at Westdahl volcano, Alaska, modeled from satellite radar interferometry[END_REF][START_REF] Lu | Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 1. Intereruption deformation, 1997-2008[END_REF][START_REF] Sturkell | Volcano geodesy and magma dynamics in Iceland[END_REF]. This behavior with exponentially decaying uplift following an eruption was monitored by interferometric synthetic aperture radar (InSAR) measurements after the 1991 Westdahl eruption [START_REF] Lu | Magma supply dynamics at Westdahl volcano, Alaska, modeled from satellite radar interferometry[END_REF] and after the 1997 eruption at Okmok Volcano [START_REF] Lu | Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 1. Intereruption deformation, 1997-2008[END_REF]. This behavior is consistent with the replenishing of a shallow reservoir from a deeper and constant pressure source through a conduit where Newtonian magma flows [Lengliné et al., 2008]. This model of a magma reservoir that is fed by a constant pressure source has also been coupled to a damage model, to interpret the seismicity rate at several basaltic volcanoes [Lengliné et al., 2008]. Time constants of one to several years have been inferred. [START_REF] Pinel | On the relationship between cycles of eruptive activity and volcanic edifice growth[END_REF] applied a similar model with a constant and deep pressure source to discuss long-term eruptive rates. However, based on ocean bottom pressure recorder measurements, [START_REF] Nooner | Volcanic inflation measured in the caldera of Axial Seamount: Implications for magma supply and future eruptions[END_REF] showed two successive trends of displacements with different time scales that followed a submarine eruption at Axial Seamount: short-term exponentially decaying uplift with a time constant in the range of 21-48 days, followed by a long-term trend, with linear or exponentially decaying inflation with a long time constant. They presented the short-term exponential decaying uplift as being specific of seamount behavior, probably as the observed time constant was smaller than that measured at surface volcanoes [Lengliné et al., 2008]. [START_REF] Nooner | Volcanic inflation measured in the caldera of Axial Seamount: Implications for magma supply and future eruptions[END_REF] proposed two potential interpretations for the short-term exponential behavior: an influx of magma from other shallow satellite magma bodies or viscoelastic relaxation plus a porous flow of melt from the partial melt region that underlies the magma chamber. It can nevertheless be argued that replenishment from a deep source with a strong hydraulic connection remains a potential explanation for this short-term behavior.

Journal of Geophysical Research: Solid Earth

10.1002/2013JB010569

Except for [START_REF] Melnik | special volume "The Eruption of Soufriere Hills Volcano, Montserrat from 2000 to 2010[END_REF], all of the studies that have dealt with temporal evolution of a displacement field associated with dynamic flow of magma into a shallow storage zone have always considered only one given reservoir that is fed by a constant pressure deep source. However, many recent deformation studies that have been based on in situ measurements at Soufriere Hills Volcano [START_REF] Elsworth | Implications of magma transfer between multiple reservoirs on eruption cycling[END_REF][START_REF] Foroozan | Dual reservoir structure at Soufrière Hills Volcano inferred from continuous GPS observations an heterogeneous elastic modeling[END_REF][START_REF] Foroozan | Magmatic metering controls the stopping and restarting of eruptions[END_REF][START_REF] Hautmann | Pressure changes in the magmatic system during the December 2008/January 2009 extrusion event at Soufrière Hills Volcano, Montserrat (W. I.), derived from strain data analysis[END_REF][START_REF] Melnik | special volume "The Eruption of Soufriere Hills Volcano, Montserrat from 2000 to 2010[END_REF] or InSAR data at Fernandina Volcano [START_REF] Chadwick | The May 2005 eruption of Fernandina volcano, Galápagos: The first circumferential dike intrusion observed by GPS and InSAR[END_REF], Galapagos Islands [START_REF] Bagnardi | Space-geodetic evidence for multiple magma reservoirs and subvolcanic lateral intrusions at Fernandina Volcano, Galápagos Islands[END_REF], have revealed that the displacement field can only be satisfactorily interpreted by considering two magma reservoirs that are located on the same magma vertical path but at two different depths. On the other hand, all of the studies that have interpreted deformation data with two magma reservoirs and with connecting conduits or dikes have inverted the displacement fields to infer each source of deformation without taking into account the existing flow dynamics between these sources.

In this study, we analyze the large ground deformation that was measured by Global Positioning System (GPS) recorders around Grímsvötn Volcano (Iceland) after the three last recent eruptions (1998, 2004, and 2011). Due to the regional geodynamic context, which is detailed below in section 2.1, this ground deformation certainly originates from tectonic, isostatic, and volcanic sources. As detailed in section 3, the first step of the study aims to estimate and define the limits of the tectonic and isostatic components of the deformation, to evaluate the remaining deformation signal that is produced by the volcanic activity of Grímsvötn (section 4). In a second step, we propose an analytical model of surface ground deformation produced during intereruptive periods by a magmatic system of two connected reservoirs (section 5). Our model considers two reservoirs at different depths on the same vertical line and takes into account the magma flow between them. In section 6, the comparison with observed GPS data allows us to add additional constraints on the reservoir size and magma inflow rates at the bottom of the volcanic system. Low-temperature conditions during a large part of the year at Grímsvötn result in icing of the GPS antenna, which creates strong seasonal perturbation to the vertical component and substantially lowers its accuracy, compared to the horizontal component. Moreover, isostatic effects also have less effect on the horizontal components, as the GPS station is located at the center of the ice cap. When comparing the GPS time series with our models, we thus only consider the horizontal components of displacement, with a limit on the vertical displacement rate used to narrow down the problem.

Geodynamic Context and Volcanic Activity at Grímsvötn

Geodynamic Context

Iceland is the emerged part of the Mid-Atlantic Ridge, which interacts with a plume at this specific location. This coupling of mid-oceanic extension and plume activity leads to intense volcanic activity that is associated with an unusual crustal thickness, which reaches 40 km at the apex of the hot spot [START_REF] Darbyshire | Structure of the crust and uppermost mantle of Iceland from a combined seismic and gravity study[END_REF]. The boundary between the North American plate and the Eurasian plate is characterized by an extension rate of 19. 4-19.6 mm/yr in a 104.5 • N direction [START_REF] Demets | Geologically current plate motions[END_REF], and they are located along the fissure swarms that form the icelandic rift zone (Figure 1).

Iceland is also covered by several ice caps, the largest of which are Vatnajökull, Hofsjökull, Langjökull, and Myrdalsjökull. Among these, with an ice volume of 3100 km 3 [START_REF] Björnsson | Icelandic glaciers[END_REF], Vatnajökull (Figure 1) is the largest ice cap in Europe, although it has been decreasing in volume rapidly since 1890, from the time of the end of the Little Ice Age. This decrease is likely to have accelerated nowadays due to global warming [START_REF] Björnsson | Icelandic glaciers[END_REF]. The overall retreat of the icelandic glaciers is causing a mechanical response in the Earth crust [START_REF] Pagli | Glacio-isostatic deformation around the Vatnajökull icecap, Iceland, induced by recent climate warming: GPS observations and finite element modeling[END_REF][START_REF] Arnadóttir | Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland[END_REF][START_REF] Pinel | Discriminating volcano deformation due to magma movements and variable surface loads: Application to Katla subglacial volcano, Iceland[END_REF]. As a first measurable consequence, an uplift of over 20 mm/yr is currently occurring around Vatnajökull ice cap [START_REF] Arnadóttir | Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland[END_REF]. The Vatnajökull ice cap covers several active volcanoes. Among these, Grímsvötn and Bárðarbunga, and their associated fissure swarms, Laki and Veidivötn, respectively, have been the most active in Iceland over the last several hundred years [START_REF] Larsen | Eight centuries of periodic volcanism at the center of the Icelandic hotspot revealed by glacier tephrostratigraphy[END_REF]].

Grímsvötn Volcano

Grímsvötn Volcano (64 1) is a subglacial basaltic volcano that is located beneath the Vatnajökull ice cap. It appears as a 10-12 km wide and 200-300 m deep, clover-shaped caldera complex. This caldera was revealed by radio-echo soundings [START_REF] Björnsson | Volcanoes beneath Vatnajökull, Iceland: Evidence from radio echo-sounding, earthquakes and jökulhlaups[END_REF], and it has been divided into three regions: the east, north, and main (or south) subcalderas [START_REF] Gudmundsson | Gravity and magnetic studies of the subglacial Grímsvötn volcano, Iceland: Implications for crustal and thermal structure[END_REF]. A shallow magma chamber has been shown by geodetic measurements [START_REF] Sturkell | Deformation of Grímsvötn volcano, Iceland: 1998 eruption and subsequent inflation[END_REF][START_REF] Sturkell | Volcano geodesy and magma dynamics in Iceland[END_REF].

• 24 ′ 29 ′′ N, 17 • 16 ′ 22 ′′ W) (Figure
Figure 1. Map of Iceland, showing the outline of the volcanic zone along the central axis at the limit of the plate (gray dashed line). The spreading rate is 9.8 mm/yr for each plate [from DeMets et al., 2010]. The Icelandic volcanic zone begins with the prolongation of Reykjanes Ridge, which is known as the Reykjanes Peninsula Rift (RPR). This extends northward by the Western Volcanic Zone (WVZ). A second fissure swarm starts with the Eastern Volcanic Zone (EVZ) and continues more to the north with the Northern Volcanic Zone (NVZ). Two transform zones are indicated (black dashed lines): the South Iceland Seismic Zone (SISZ), which links the WVZ and EVZ, and the Tjornes Fracture Zone (TFZ) which links the NVZ and Kolbeinsey Ridge. The dashed box shows the study area, with Grímsvötn Volcano indicated (G) (64 • 24 ′ 29 ′′ N, 17 • 16 ′ 22 ′′ W), under the Vatnajökull ice cap (8000 km 2 ), with the position of Laki fissure also indicated. From [START_REF] Sturkell | Volcano geodesy and magma dynamics in Iceland[END_REF].

Also, based on tomographic inversion of micro-earthquake activity, [START_REF] Alfaro | Structure of the Grímsvötn central volcano under the Vatnajökull icecap, Iceland[END_REF] estimated that this reservoir is a sill that is located ∼3 km to 4 km in depth, with a maximum extent of 7 km to 8 km E-W and 4 km to 5 km N-S, and with a thickness of 1 km. From the recent inversion of the ground deformation data (using GPS and a tiltmeter), the top of the chamber is expected to be located at an even shallower depth, at ∼1.7 km [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF].

Grímsvötn Volcano is the most active volcano in Iceland, with an eruption history that has been characterized by distinct 50 year long to 80 year long periods of high eruption rates alternating with equally long periods of low eruption rates [START_REF] Larsen | Eight centuries of periodic volcanism at the center of the Icelandic hotspot revealed by glacier tephrostratigraphy[END_REF]. For the last 200 years (1816-2011), the eruptive history has been well documented, with 21 events recorded [START_REF] Thordarson | Volcanism in Iceland in historical time: Volcano types, eruption styles and eruptive history[END_REF]. A new cycle of high activity began in 1996, with a flank eruption in Gjálp, and then three eruptions occurred, one each in [START_REF] Larsen | Eight centuries of periodic volcanism at the center of the Icelandic hotspot revealed by glacier tephrostratigraphy[END_REF][START_REF] Zeeuw-Van Dalfsen | Satellite radar interferometry 1993-1999 suggests deep accumulation of magma near the crust-mantle boundary at the Krafla volcanic system, Iceland[END_REF]and 2011[START_REF] Oddsson | Monitoring of the plume from the basaltic phreatomagmatic 2004 Grímsvötn eruption-application of weather radar and comparison with plume models[END_REF] estimated a volume of erupted magma of 0.02 ± 0.004 km 3 dense rock equivalents (DRE) for the 2004 eruption, and [START_REF] Gudmundsson | The May 2011 eruption of Grímsvötn[END_REF] calculated a volume of 0.2 km 3 to 0.3 km 3 DRE for the 2011 eruption.

GPS Analysis and Data Processing

GPS Data Analysis

Daily site positions were obtained using the International Global Navigation Satellite System Service (IGS) final precise orbits [START_REF] Beutler | The International GPS Service (IGS): An interdisciplinary service in support of Earth sciences[END_REF], as well as IGS Earth rotation parameters and data from nearby permanent GPS stations. Absolute antenna phase center offset models were used during this analysis [START_REF] Bilich | Three-method absolute antenna calibration comparison[END_REF]. The data from the permanent GPS network for the 2000 to 2012 period (Icelandic GPS stations: DYNC, EYVI, FJOC, HAUC, KIDC, KVSK, GFUM, SKRO, and STKA; and IGS stations: HOFN, KELY, KIR0, KIRU, KOSG, KULU, MAR6, NAIN, NYA1, ONSA, QAQ1, THU1, TRO1, REYK, STJO, and WSRT) were analyzed by Bernese 5.2 software using the following strategy: (1) initial ionosphere-free analysis with computation of residuals; (2) residuals analysis; (3) resolution of ambiguities using the quasi-ionosphere-free strategy of resolution, with the ionosphere model obtained from Center Orbit Determination in Europe; and (4) computation of the normal equations. The troposphere-induced propagation delays were estimated from observations made every 2 h, using the Saastamoinen model [START_REF] Saastamoinen | Contribution to the theory of atmospheric refraction[END_REF] to correct for the atmospheric delays, and the NIELL mapping function. Each daily solution is transformed into the International Terrestrial Reference Frame (ITRF) 2008 [START_REF] Altamimi | ITRF2008 plate motion model[END_REF] with a six-parameter Helmert solution (three translation parameters, three rotation parameters) using the IGS stations: KELY, KIR0, KIRU, KOSG, KULU, MAR6, NAIN, NYA1, ONSA, QAQ1, THU1, TRO1, STJO, and WSRT. Discontinuities originally in the ITRF2008 solutions are introduced. Outlier detection is performed using the Bernese 5.0 software, with the ITRF2008 site-weighted coordinates, and the velocities constrained at their ITRF2008 values for points defined in ITRF2008.

Local effects, such as ice-covered GPS antennae, episodic loading by snow, or rainfall, clearly affect the recorded signals, and particularly the vertical component, which often shows seasonal effects [START_REF] Geirsson | Current plate movements across the Mid-Atlantic Ridge determined from 5 years of continuous GPS measurements in Iceland[END_REF][START_REF] Grapenthin | Icelandic rhythmics: Annual modulation of land elevation and plate spreading by snow load[END_REF]. Velocities and time series expressed in the ITRF2008 reference were finally expressed in the reference framework of the plate limit of the North America-Eurasia plates, defined using the ITRF2005/North America and ITRF2005/Eurasia rotation poles [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF]. The 10 stations are arranged several tens of kilometers from each other (see Figure 2 for station positions). These GPS stations were established at different times, with some recording from as early as 1997 (e.g., HOFN station) and others from 2008 (e.g., HAUC station). From a temporal point of view, the coverage is relatively good for these stations (see Figures S1 andS2 in the supporting information) for the time interval between 2006 and 2012.

Various Contributions to the Ground Deformation

The study area (Figure 1) might encompass different sources of displacement of tectonic, isostatic, volcanic, or local origins. The tectonic contribution has been estimated at ∼1 cm/yr in previous studies [START_REF] Demets | Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions[END_REF][START_REF] Arnadóttir | Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland[END_REF]. With this source of displacement characterized by a long time scale, this tectonic effect can be considered as steady over a 5 year study, and this is assimilated into a linear trend in our records. Vatnajökull ice cap has been retreating since 1890, as have all of the Icelandic ice caps, which has induced an almost linear vertical uplift of ∼20 mm/yr around the ice cap that is associated with a horizontal and divergent displacement field of a smaller amplitude (≃3-4 mm/yr [START_REF] Pagli | Glacio-isostatic deformation around the Vatnajökull icecap, Iceland, induced by recent climate warming: GPS observations and finite element modeling[END_REF]). Due to the presence of active volcanoes, such as Grímsvötn and Bárðarbunga, volcanic sources also have major roles in the ground deformation [START_REF] Sturkell | Volcano geodesy and magma dynamics in Iceland[END_REF].

The Tectonic Component of the Displacement and the Isostatic Rebound Effect

We analyzed here the velocities recorded at the GPS stations, except at station GFUM, which was clearly affected by the eruptive activity. For the vertical component of all the stations (Figure 2B), the displacement is upward and shows a linear trend (supporting information Figures S1 andS2) at a large rate (≃20 mm/yr), and this can be attributed to glacial isostatic rebound [START_REF] Pagli | Glacio-isostatic deformation around the Vatnajökull icecap, Iceland, induced by recent climate warming: GPS observations and finite element modeling[END_REF][START_REF] Auriac | Iceland rising: Solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling[END_REF][START_REF] Arnadóttir | Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland[END_REF]. For the horizontal components, for most of the stations except DYNC, the displacement also follows Journal of Geophysical Research: Solid Earth 10.1002/2013JB010569 Figure 3. Observed and modeled horizontal component displacements. The plate limit is approximated by a vertical dike of infinite length and locked at 8 km in depth, in an elastic half-space. The observed horizontal displacements are corrected for a radial GIA of 5 mm/yr amplitude at all of the stations except GFUM with a center in the middle of the ice cap (point V C ). The best model is obtained for an opening of 20 mm/yr with 8 mm/yr of sinistral strike-slip component.

a clear linear trend. The magnitudes of the displacement rate are almost 10 mm/yr for all of the western stations and larger for the eastern stations, at about 15 mm/yr. Furthermore, the directions of the velocities are, as expected, almost perpendicular to the rift zone and are opposite for the GPS stations located on the Eurasian plate compared to those located on the North American plate. [START_REF] Demets | Geologically current plate motions[END_REF] estimated the half-spreading rate at 9.7 mm/yr to 9.8 mm/yr, and [START_REF] Arnadóttir | Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland[END_REF] simulated a rate of full opening for EVZ at 21 ±2 mm/yr. From these observations, we can attribute the observed linear trend of the time series to the tectonic component and the Glacial Isostatic Adjustment (GIA).

Volcanic Deformation at GFUM

GFUM station is located 3.5 ± 0.2 km from the center of Grímsvötn caldera [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF], and it is fixed on the Nunatak Svíahnukur Eystri. At GFUM (Figure 2), the horizontal displacement rate (≃50 mm/yr) is fivefold greater than the rate observed at the other stations, which clearly indicates a major volcanic component in the observed ground deformation. To separate the volcanic and tectonic components, we estimated the tectonic component at GFUM by modeling the rift opening rate and direction with a dislocation model [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. In this model, we simulated the EVZ by an infinite and vertical dike, with a direction 43 • N (Figure 3) in a uniform elastic half-space. During the inversion, all of the continuous GPS velocities, except that of GFUM, were considered. To be able to compare horizontal velocities free of isostatic rebound, we decided to estimate and subtract a radial divergent displacement from our data set, with a center in the middle of Vatnajökull ice cap. We selected the best GIA correction as the correction that allows the best numerical simulation for the tectonic component (see supporting information Table S1). The best solution is obtained for a 5 ± 1 mm/yr GIA, which is slightly higher than the estimation of [START_REF] Pagli | Glacio-isostatic deformation around the Vatnajökull icecap, Iceland, induced by recent climate warming: GPS observations and finite element modeling[END_REF]. We thus subtracted 5 mm/yr from the observed linear displacement trend. The corrected linear trends are reported in Figure 3. The best solution is obtained with an oblique opening rate: 20 mm/yr for the normal component and 8 mm/yr for the sinistral strike-slip component, locked at 8 km in depth (see Figure 3). It must be underlined that the corrected DYNC time series presents a null horizontal velocity, which is consistent with its location on the limit of the plate. We predict a tectonic velocity for GFUM with a 7.5 mm/yr east and -2.7 mm/yr north component relative to ridge-fixed reference frame. Horizontal components of the GFUM time series are then corrected for the tectonic contribution by removing 8.0 mm/yr from the slope of the time series (Table 1). The horizontal component of the GIA at GFUM can be neglected because it is close to 0 in the middle of the ice cap. Hereafter, we consider 42 mm/yr for the horizontal volcanic residual trend at GFUM,

The Volcanic Component of the GPS Displacement at GFUM Station

Figure 4 shows the displacement for each of the components at GFUM station during the complete intereruptive cycle between the 2004 and 2011 eruptions. Figure 5a shows the direction of the coeruptive and posteruptive displacement from this study for the 2004 and 2011 eruptions and from [START_REF] Sturkell | Deformation of Grímsvötn volcano, Iceland: 1998 eruption and subsequent inflation[END_REF] for the 1998 eruption (seven field campaign measurements made at GRIM, located 250 m away from the continuous GPS station GFUM). We note almost the same N-W direction (330 • N for 2004 eruption and 320 • N for 2011 eruption) for the syneruptive displacements and a strictly opposite direction (150 • N and 140 • N for 2004 and 2011 eruption, respectively) for the posteruption displacements. The important feature to note is the stability of the direction (almost 145 • N) of the displacement for each of the three eruptions, despite the variation in the event location, which indicates that the source of deformation remained at the same location during these last three eruptions. During the 2004 eruption, the observed coeruptive displacement reached 11 ±1 cm for the east component, 24 ±1 cm for the north component, and 20 ±1 cm for the vertical component (Figure 4). After the 2004 eruption, a fast readjustment occurred that lasted for a few months and was followed by a linear trend in the three components. The displacement for the 325 • N direction, called the radial displacement in the following, is calculated and represented in Figure 5b. After the 2011 eruption, similar behavior is observed for the radial displacement as following the previous eruption. This result is confirmed in Figure 5b, where the radial displacement observed at the Grímsvötn GPS station during the last three posteruptive periods shows strikingly similar features: an initial exponential behavior, followed by linear evolution. This displacement is consistent with magma replenishment within the shallow reservoir modeled by [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF]; however, it is controlled by two processes.

Journal of Geophysical

A Two-Magma-Chamber Model

In this section, we describe an analytical model with two connected magmatic reservoirs that are filled at a constant basal rate (see Figure 6), and we express the pressure changes in the two reservoirs, the cumulative volume stored at a shallow level, and the ground deformation produced at the surface.

Analytical Model

The displacement observed at the GFUM GPS station during the last three coeruptive and posteruptive periods shows constant direction for the movement and the same temporal behavior. For the temporal behavior, the three posteruption inflation stages can all be decomposed into an exponential phase of similar duration, of ∼1 year, followed by a linear trend (Figure 5b). This behavior cannot be explained by a model of a single shallow magmatic source filled from an infinite deep magmatic reservoir at constant pressure. Indeed, such a model would give an exponential trend tending asymptotically to a constant value when the pressure difference between the two reservoirs vanishes [see e.g., Lengliné et al., 2008;[START_REF] Segall | Volcano deformation and eruption forecasting[END_REF]. In this model, the magma is characterized by its viscosity 𝜇 and density 𝜌 m and is supplied as a constant flux Q in to the deep reservoir. The magmatic flux between the two reservoirs is denoted Q. 𝜌 r is the density of the surrounding rock. P s (t) and P d (t) are the shallow and deep reservoir pressures, respectively, which evolve through time.

Journal of Geophysical

To explain the linear trend following the exponential trend observed at Grímsvötn, we propose a model with two finite reservoirs in which the deep reservoir is continuously filled at the bottom at a constant rate Q in . The geometry of this model is presented in Figure 6.

We consider two possible geometries, as spherical or sill-like, for the two magmatic reservoirs that are located below Grímsvötn Volcano and are embedded in an infinite elastic half-space. They are connected by an open conduit of radius a c and length H c . We assume that the ascent of magma between the two reservoirs can be described as a Poiseuille flow of rate Q [START_REF] Pinel | Magma chamber behavior beneath a volcanic edifice[END_REF].

Q = 𝜋a c 4 8𝜇 [ - dP dz -𝜌 m g ] , (1)
where dP dz is the vertical hydrostatic pressure gradient, g is gravity, z is the vertical coordinate pointing upward, and 𝜇 and 𝜌 m are the viscosity and density of the magma, respectively, which are assumed constant throughout the system. The magma is thus considered as a first approximation as an incompressible fluid, such that in the following, the mass balance is directly derived from the volume balance.

Expression of Overpressure and Magma Inflow for the Two Reservoirs

The pressure in the deep reservoir P d (t) is the sum of the lithostatic pressure at the reservoir depth H d and an overpressure ΔP d (t):

P d (t) = P litho (-H d ) + ΔP d (t) = 𝜌 r gH d + ΔP d (t), (2)
where 𝜌 r is the density of the rock in the medium. Similarly, the pressure in the shallow reservoir, at depth H s , can be written as the sum of the lithostatic pressure and an overpressure ΔP s (t):

P s (t) = P litho (-H s ) + ΔP s (t)
= 𝜌 r gH s + ΔP s (t).

(3)

From equations ( 1)-(3) for the flow Q from the deep to the shallow reservoir, we obtain

Q(t) = 𝜋a c 4 8𝜇 [ P d (t) -P s (t) H c -𝜌 m g ] = 𝜋a c 4 8𝜇H c [(𝜌 r -𝜌 m )gH c + ΔP d (t) -ΔP s (t)]. (4)
Journal of Geophysical Research: Solid Earth

10.1002/2013JB010569

As the deep reservoir is fed with magma at a flow rate of Q in and transfers magma to the shallow reservoir at a flow rate of Q, the mass balance for the deep reservoir can be written as

dΔV d (t) dt = Q in -Q(t), ( 5
)
where ΔV d is the volume of magma injected into the deep reservoir. The mass balance for the shallow reservoir can likewise be written as

dΔV s (t) dt = Q(t), ( 6
)
where ΔV s is the volume of magma injected into the shallow reservoir. After [START_REF] Delaney | Volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kilauea volcano[END_REF], the relations between the injected volume and the overpressure variation in the two reservoirs are

ΔV d (t) = ΔP d (t) 𝛾 d 𝜋a d 3 G (7)
for the deep reservoir and

ΔV s (t) = ΔP s (t) 𝛾 s 𝜋a s 3 G (8)
for the shallow one, where G is the rigidity modulus, and

𝛾 = 8 3𝜋 (1 -𝜈) (9)
for a sill (oblate reservoir with a small thickness with regard to its lateral dimension) [START_REF] Amoruso | Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data[END_REF] or 𝛾 = 1 for a spherical reservoir [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them[END_REF].

As detailed in Appendix A, from equations ( 4), ( 6), and ( 8), we obtain the first-order linear differential equation for ΔP s (t):

dΔP s (t) dt + 𝜉ΔP s (t) = Ga c 4 8𝛾 s 𝜇H c a s 3 [ ΔP d (0) + 𝛾 s k𝛾 d ΔP s (0) + (𝜌 r -𝜌 m )gH c + G 𝛾 d 𝜋a d 3 Q in t ] (10)
with

𝜉 = Ga c 4 (𝛾 s + 𝛾 d k) 8𝜇H c a s 3 𝛾 s 𝛾 d k ,
and k is the ratio of the reservoir radius to the power of 3:

k = a d 3 a s 3 .
The solution to this differential equation can be expressed as

ΔP s (t) = A(1 -e -𝜉t ) + GQ in 𝜋a s 3 (𝛾 s + 𝛾 d k) t + ΔP s (0) (11)
with

A = 𝛾 d k 𝛾 s + 𝛾 d k [ ΔP d (0) -ΔP s (0) + (𝜌 r -𝜌 m )gH c - 8𝛾 s Q in 𝜇H c 𝜋a c 4 (𝛾 s + 𝛾 d k)
] .

Substituting ΔP s (t) from equation (11) into equation (A6), we get the overpressure in the deep reservoir ΔP d (t) as a function of time:

ΔP d (t) = - 𝛾 s A 𝛾 d k (1 -e -𝜉t ) + GQ in 𝜋a s 3 (𝛾 s + 𝛾 d k) t + ΔP d (0). ( 12
)
Both pressure functions ΔP s (t) and ΔP d (t) are the sum of an exponential function and a linear function.

Their evolution is illustrated in Figure 7. Just after an eruption, the overpressure in the shallow reservoir rises exponentially, whereas the overpressure in the deep reservoir decreases exponentially, both with the same time constant 1∕𝜉. Then, both pressures increase with the same linear trend. The slope of the linear trend depends on the basal magma inflow, the crustal rigidity, and the geometry of both of the magma reservoirs. in the shallow and deep reservoirs, respectively. After an eruption, the overpressure within the shallow reservoir increases exponentially, while in the deep reservoir, the overpressure decreases (this decrease depends on the ratio a d ∕a s between the reservoir sizes). Then, after a given time that is characterized by a time constant 𝜒, the two overpressures increase at the same rate. The time origin (t = 0) is set just after the eruption.

Journal of Geophysical

After equation ( 12), the amplitude of the exponential decrease of the deep reservoir overpressure ΔP d (t) depends on k, i.e., on the relative sizes of the reservoirs. The larger the deep reservoir is, the lower its posteruption pressure drop. The differential pressure between the two reservoirs tends to a constant when the exponential terms become negligible.

As described before, the magma input rate Q in into the system is chosen as constant, to match the linear trend of the ground surface deformation observed at Grímsvötn. The magma transfer rate Q(t) from the deep to the shallow reservoir is obtained from equations ( 11), ( 12), and ( 4) :

Q(t) = A 𝜋a 4 c 8𝜇H c 𝛾 s + 𝛾 d k 𝛾 d k e -𝜉t + 𝛾 s Q in 𝛾 s + 𝛾 d k . ( 13
)
When t tends to infinity,

Q(t) tends to Q ∞ : Q ∞ = 𝛾 s Q in 𝛾 s + 𝛾 d k . ( 14
)
Figure 8 shows Q ∞ ∕Q in , which is the ratio between the volume rate transferred to the shallow reservoir and the volume rate entering into the deep reservoir. For a large radius of the deep reservoir compared to the shallow reservoir, the magma is mainly stored in the deep reservoir, and a small proportion of the input magma is transferred to the surface. The cumulative volume that enters into the shallow reservoir between time 0 and t is obtained by the integration of equation ( 13):

ΔV(t) = Q ∞ t + A 𝜋a 3 s 𝛾 s G (1 -e -𝜉t ). ( 15
)
Figure 8. Evolution of the flux ratio Q ∞ ∕Q in as a function of the radius a d ∕a s for a shallow sill-like reservoir and a sill-like or spherical deep reservoir. Q in is the basal flux of magma that enters the deep reservoir (see Figure 6), and Q ∞ is the steady magma flux transferred to the shallow reservoir after the reequilibration stage.

At an observation time t f which is large enough compared to the exponential time constant (for t f >> (1∕𝜉)), this expression can be approximated by

ΔV(t f ) =ΔV cont + ΔV temp ( 16
)
=Q ∞ t f + A 𝜋a 3 s 𝛾 s G .
The term ΔV temp = A(𝜋a 3 s 𝛾 s )∕G corresponds to the volume due to the replenishment of the shallow reservoir by the drawing up of magma from the deep source after the eruption, in a posteruptive reequilibration stage, whereas Δ V cont = Q ∞ t f corresponds to the volume due to continuous basal filling at a rate Q in .

Expression of Ground Deformation

The total displacement vector at the surface can be obtained by adding together the displacement vectors created by the deep and shallow pressure sources. In the following, we consider the cases in which each reservoir can be either a sphere or a sill.

Journal of Geophysical Research: Solid Earth

10.1002/2013JB010569

For a reservoir with an overpressure ΔP, the horizontal, radial to the source U R (t), and vertical U V (t) surface displacement, was given by [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them[END_REF] for a sphere and is derived for a sill (Appendix B):

U R (t) = (1 -𝜈)r G 𝛼 ( a 3 R 3 ) ΔP(t) (17) U V (t) = U R (t) H r , ( 18
)
where R is the distance from the reservoir to the surface station (Figure 6), r is the radial distance between the source axes and the GPS station, and 𝜈 is the Poisson coefficient. For a spherical reservoir 𝛼 = 1, and 𝛼 = 4H 2 ∕𝜋R 2 for a sill-like reservoir, where H is the depth of the reservoir.

The horizontal surface displacement produced at the surface by a shallow reservoir with overpressure ΔP s (t) and a deep reservoir with overpressure ΔP d (t) is given by

U R (t) = (1 -𝜈) G r ( 𝛼 s a 3 s R s 3 ΔP s (t) + 𝛼 d a 3 d R d 3 ΔP d (t) ) , ( 19
)
where R s and R d are the distances from the shallow and the deep reservoir to the surface station, respectively, and H s and H d are the depths of the shallow and deep reservoirs, respectively (Figure 6).

Substituting the expressions of ΔP s (t) and ΔP d (t) from equations ( 11) and ( 12) into equation ( 19) gives

U R (t) = A(1 -e -𝜉t ) (1 -𝜈)r G ( 𝛼 s a 3 s R s 3 -𝛼 d 𝛾 s a 3 d k𝛾 d R d 3 ) + UR∞ t, ( 20
)
with UR∞ , the linear horizontal trend, as

UR∞ = (1 -𝜈) 𝜋a 3 s 𝛾 s r ( 𝛼 s a 3 s R s 3 + 𝛼 d a 3 d R d 3 ) Q ∞ . ( 21
)
We can estimate the vertical displacement rate UV∞ as

UV∞ = (1 -𝜈) 𝜋a 3 s 𝛾 s ( 𝛼 s H s a 3 s R s 3 + 𝛼 d H d a 3 d R d 3 ) Q ∞ . ( 22
)
At an observation time t f that is large enough compared to the exponential time constant (i.e., for t f >> (1∕𝜉)), the horizontal surface displacement can be expressed as a function of the two contributive volumes ΔV cont and ΔV temp :

U R (t f ) = (1 -𝜈)r𝛼 s 𝜋R s 3 𝛾 s [ ΔV temp ( 1 -𝜖 𝛾 s k𝛾 d ) + ΔV cont (1 + 𝜖) ] , ( 23
)
with

𝜖 = 𝛼 d 𝛼 s R 3 s R 3 d a 3 d a 3 s .
In our particular case, the lateral distance of the GPS station is r = 3.5 km, according to the location of GFUM station with regard to the caldera center (see Figure 5a). As 𝜖∕k = R 3 s ∕R 3 d , and H s ≤4 km [START_REF] Alfaro | Structure of the Grímsvötn central volcano under the Vatnajökull icecap, Iceland[END_REF], we can consider that for the deep reservoir at a depth >10 km, the term 𝜖 𝛾 s k𝛾 d can be neglected, such that equation ( 23) becomes

U R (t f ) = (1 -𝜈)r𝛼 s 𝜋R s 3 𝛾 s [ΔV temp + ΔV cont (1 + 𝜖)]. ( 24
)
Using equations ( 15) and ( 23), we can deduce that the cumulated volume of input magma into the shallow reservoir at a time t f is given by

ΔV(t f ) = UR∞ t f 𝜋R s 3 𝛾 s (1 -𝜈)r𝛼 s U R (t f ) UR∞ t f - 𝜖 1 + 𝜖 ] . ( 25
)
Journal of Geophysical Research: Solid Earth 10.1002/2013JB010569

Estimation of the Model Parameters

In this section, we aim to estimate the model parameters through comparison of the horizontal and vertical displacement observed at Grímsvötn with the analytical model. Because of the uncertain GIA contribution, the intense icing perturbation at GFUM, and the lower accuracy of the vertical component compared to horizontal components, we use the horizontal displacements time series at GFUM to constrain the parameters of the model and we use the cumulative vertical displacement as a limit for the solutions (Figure 4).

We first focus on the exponential rise in the displacement following the eruptions. From equation ( 17), the displacement U R (t) is linearly related to the overpressure ΔP s (t) and ΔP d (t) in the reservoirs. The exponential time constant 𝜏 = 1∕𝜉 in ΔP s (t) and ΔP d (t) can thus be inferred from U R (t). Figure 5b shows the horizontal displacement observed at GFUM during the last three posteruptive periods, matched by a function of the type Φ(1 -e -t∕𝜏 ) + U∞ t + C (20). The best exponential fits give 𝜏 = 0.33 ± 0.08 year and 𝜏 = 0.13 ± 0.04 year for the 2004 and 2011 posteruptive periods, respectively. Due to the small number of measurements available for the 1998 posteruptive period, the time constant was not estimated. The two time constants estimated differ by 2 months, which can be explained by a small change in the conduit geometry or effects of magma compressibility. However, these time constants are significantly smaller than those obtained for other basaltic volcanoes. Lengliné et al. [2008] reported 𝜏 = 1.13 years for Kilauea Volcano and 𝜏 = 7.3 years for Mauna Loa Volcano, and [START_REF] Lu | Magma supply dynamics at Westdahl volcano, Alaska, modeled from satellite radar interferometry[END_REF] reported 𝜏 = 6 years for Westdahl Volcano. This difference can easily be explained either by a stronger hydraulic connection between the two reservoirs (a larger conduit with a smaller vertical extension) or by less viscous magma. The magma viscosity at Grímsvötn has been measured as between 8 Pa s and 2000 Pa s [START_REF] Hobiger | Viscosity characteristics of selected volcanic rock melts[END_REF], which is in the same range as other basaltic volcanoes that have been studied; thus, a different conduit geometry is expected. It is of note that the time constant measured at Grímsvötn is close to that for the short-term response at the submarine Axial Seamount Volcano given by Nooner and Chadwick [2009], which indicates a common, more efficient magma transfer in a ridge context. Now we focus on the linear trend of the displacement. This linear trend UR∞ that is observed on the U R displacement when t → ∞ during the post 2004 eruption period is the best documented. During the post-2011 eruption period, the trend is almost linear but with a larger slope (60 mm/yr) than for the post-2004 one (42 mm/yr). This could be due to a larger basal magma inflow Q in , after a 2011 eruption where the erupted magma volume was 10 times larger than for the 2004 eruption. We used the value obtained for the post-2004 period, i.e., UR∞ for the following. For a given geometry model, we used the values of the linear trend UR∞ and total horizontal displacement U R (t f ) to estimate the vertical displacement UV ∞ and the cumulative volume ΔV, according to equations ( 22) and ( 25). From equations ( 21) and ( 14), we can also estimate the volume rates Q ∞ and Q in .

Figure 9 shows the results of the modeling. The geometry of the shallow reservoir is set as a sill, after [START_REF] Alfaro | Structure of the Grímsvötn central volcano under the Vatnajökull icecap, Iceland[END_REF]. Its depth was inferred at 1.7 km (±0.2 km) by inversion of the horizontal and vertical components of the GPS station together with the tiltmeter data using a Mogi model [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF]. Using the same inversion for a sill, the inferred depth should be multiplied by 5 3 (combining equations derived in [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF] and equations (B4) and ( B5)), such that we fix the shallow reservoir depth H s at 3 km. For the deep reservoir, we consider two possible geometries: sphere or sill. The two parameters of the model are the depth of the deep reservoir H d , which varies from 10 km to 35 km, and the reservoir radius ratio a d ∕a s , which varies from 1 to 3. For each set of parameters, we calculate the vertical displacement rate UV ∞ , the magma input rate into the deep reservoir Q in , and the cumulative volume ΔV in the shallow reservoir for t f = 6.5 years. The model validity domain (Figure 9, outside the hatched area) corresponds to values of UV ∞ <50 mm/yr, the observed rate at GFUM, which represents an upper limit without taking into account the glacial readjustment.

As UV ∞ , ΔV, and Q ∞ are linearly related to UR∞ (equations ( 21), ( 22), and ( 25)), an error in the estimation of the linear trend of horizontal displacement UR∞ will lead to the same relative error for these parameters. For instance, a change from 42 mm/yr to 43 mm/yr in UR∞ , i.e., a relative error of 5%, would give rise to the same 5% error in 6). U V ∞ is the vertical displacement rate, Q in is the basal magma flux, ΔV is the total cumulative volume in the shallow reservoir (see Figure 6 and equation ( 25)). The shallow reservoir depth (H s ) is fixed at 3 km. Nonvalid solutions (hatched area) correspond to vertical uplift rates at GFUM that are larger than 50 mm/yr. The white line corresponds to a vertical uplift rate of 50 mm/year at GFUM.

values of the parameters presented in Figure 9 do not directly depend on the shallow reservoir radius a s , but on the radii ratio a d /a s . The total volume of magma involved in the 2004 eruption was estimated to be 0.02 ± 0.004 km 3 DRE [START_REF] Oddsson | The Grímsvötn eruption in 2004: Dispersal and total mass of tephra and comparison with plume transport models[END_REF], and the volume erupted in 2011 was tenfold greater [START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF]. The larger volume of the 2011 eruption can be explained by a different magma compressibility and by a larger pressure drop, as shown by the larger coeruptive displacement. For H d >15 km, the lowest values of ΔV(t f ) range around 0.06 km 3 , which is about triple the estimation from [START_REF] Oddsson | The Grímsvötn eruption in 2004: Dispersal and total mass of tephra and comparison with plume transport models[END_REF]. The relative error on ΔV(t f ) is ∼ 70%, mainly due to the error in the location of the shallow reservoir (H s and r). With such an error, the estimated erupted volume is in the range of the model and is consistent with the fact that the volume of magma that entered the shallow reservoir after the 2004 eruption is larger than the volume emitted during this eruption, as shown by the larger amplitude of the 2004 posteruptive displacement than the 2004 coeruptive displacement. To keep only model solutions with a vertical displacement rate <50 mm/yr, the size ratio has to remain <2.5, as shown on Figure 9. As the radius of the shallow reservoir is <4 km [START_REF] Alfaro | Structure of the Grímsvötn central volcano under the Vatnajökull icecap, Iceland[END_REF], it follows that the radius of the deep reservoir has to remain <10 km. The deep reservoir is therefore focused and cannot be considered as the diffuse part of the upper mantle. Such a deep reservoir is large enough to store the 15 km 3 of the Laki eruption (1783-1785) when a vertical extension of 50 m is considered. The geometry (sphere versus sill) of the deep reservoir mainly influences the magma flux feeding the deep reservoir, which ranges between 0.01 and 0.05 km 3 /yr in the case of a sill-like deep reservoir and between 0.01 and 0.11 km 3 /yr for a spherical reservoir. If we consider, as did [START_REF] Bindeman | Time constraints on the origin of large volume basalts derived from O-isotope and trace element mineral zoning and U-series disequilibria in the Laki and Grímsvötn volcanic system[END_REF], that the same magma chamber fed the 1783-1784 AD fissure eruption of Laki (Iceland) and the Grímsvötn eruptions, and that a total volume of 15 km 3 erupted at Laki, these fluxes indicate that between 1500 years and 130 years are needed to fill the deep magma chamber. These times are clearly consistent with the residence time of the Laki magma, which was estimated at several hundreds of years from the geochemical data [START_REF] Bindeman | Time constraints on the origin of large volume basalts derived from O-isotope and trace element mineral zoning and U-series disequilibria in the Laki and Grímsvötn volcanic system[END_REF]. This several hundreds of years residence time favors an oblate geometry, which takes longer to fill than a spherical geometry. These values appear consistent with the expected magma production, which is estimated to be ∼0.17 km 3 /yr for the rift [START_REF] Pagli | Will present day glacier retreat increase volcanic activity[END_REF][START_REF] Sigmundsson | Multiple effects of surface load changes and associated stress change on magmatic systems[END_REF], and they are only slightly greater than the magma input estimated beneath Hekla Volcano, based on GPS measurements [START_REF] Geirsson | Volcano deformation at active plate boundaries: Deep magma accumulation at Hekla volcano and plate boundary deformation in south Iceland[END_REF].

Journal of Geophysical Research: Solid Earth 10.1002/2013JB010569

Discussion

The temporal evolution of surface displacement observed in volcanic areas is the result of a convolution between a history of a pressure source and the rheology of the crust and the magma. It is difficult to distinguish between these two components, such that most of the models that are used to interpret deformation measurements can be classified into two end-members. The first set of models considers a magma reservoir surrounded by an elastic homogeneous medium and fed or drawn out through a hydraulic connection [START_REF] Lu | Magma supply dynamics at Westdahl volcano, Alaska, modeled from satellite radar interferometry[END_REF][START_REF] Lu | Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 1. Intereruption deformation, 1997-2008[END_REF][START_REF] Mastin | Constraints on the size, overpressure, and volatile content of the Mount St. Helens magma system from geodetic and dome-growth measurements during the 2004-2006+ eruption[END_REF], thus favoring a complex pressure evolution within the magma reservoir. In contrast, the second set of models usually considers a more simplistic pressure scenario, most often as a sudden pressurization that acts over a viscoelastic or poroelastic medium. For instance, a sudden pressure increase in a reservoir surrounded by a viscoelastic shell in an elastic medium induces instantaneous displacement, followed by continuous displacement with a decreasing rate, following an exponential law [START_REF] Dragoni | Displacement and stress produced by a pressurized, spherical magma chamber, surrounded by a viscoelastic shell[END_REF]]. As explained by [START_REF] Dzurisin | Continuing inflation at Three Sisters volcanic center, central Oregon Cascade Range, USA, from GPS, leveling, and InSAR observations[END_REF], a shallow reservoir embedded in an elastic medium and fed by a constant pressure source through a hydraulic connection, and a sudden pressure increase in a reservoir surrounded by a viscous shell can both equally well explain exponential decay of the inflation rate. Also, thermo-poro-elastic models have been developed to understand hydrothermal systems, and they have shown that fluid injection at a constant rate can also induce surface displacement that is characterized by short-term exponential behavior followed by a linear trend [START_REF] Fournier | Understanding volcano hydrothermal unrest from geodetic observations: Insights from numerical modeling and application to White Island volcano, New Zealand[END_REF]. In this study, we have favored a model based on the simplest rheology, the elastic behavior, and the hydraulic connection between a shallow magma reservoir and a deep magma reservoir. This provides a simple and consistent explanation for the temporal evolution of the surface displacement recorded above Grímsvötn volcano.

The presence of a shallow reservoir below Grímsvötn has been revealed through several geophysical measurements [START_REF] Alfaro | Structure of the Grímsvötn central volcano under the Vatnajökull icecap, Iceland[END_REF][START_REF] Sturkell | Volcano geodesy and magma dynamics in Iceland[END_REF][START_REF] Hreinsdottir | Volcanic plume height correlated with magma-pressure change at Grímsvötn Volcano, Iceland[END_REF], which have enabled us to constrain its geometry and depth. Regarding the deep plumbing system, other volcanoes in Iceland have a deep magma chamber, such as Hekla, at about 14 km to 20 km [START_REF] Ofeigsson | Deep magma storage at Hekla volcano, Iceland, revealed by InSAR time series analysis[END_REF], or have two magma chambers, with the deep one at 21 km for Krafla [START_REF] Zeeuw-Van Dalfsen | Satellite radar interferometry 1993-1999 suggests deep accumulation of magma near the crust-mantle boundary at the Krafla volcanic system, Iceland[END_REF]. At Laki, Gjálp, and Grímsvötn, the presence of a unique same deep magma reservoir is supported by several petrological and geochemical observations [START_REF] Sigmarsson | The 1996 and 1998 subglacial eruptions beneath the Vatnajökull ice sheet in Iceland: Contrasting geochemical and geophysical inferences on magma migration[END_REF][START_REF] Steinthorsson | Petrochemistry of the Gjálp-1996 subglacial eruption, Vatnajökull, SE Iceland[END_REF][START_REF] Bindeman | Time constraints on the origin of large volume basalts derived from O-isotope and trace element mineral zoning and U-series disequilibria in the Laki and Grímsvötn volcanic system[END_REF], with, as main argument, the presence of a uniform tholeiitic major element composition of erupted basalt. Furthermore, a connection between Grímsvötn and Laki is suggested by the simultaneous activity that occurred at these two volcanoes in 1783-1784 [START_REF] Thordarson | The Laki (Skaftár Fires) and Grímsvötn eruptions in 1783-1785[END_REF]. At that time, the Laki fissure eruption produced about 15 km 3 of lava. Moreover, the uniform geochemical compositions of basalts from Grímsvötn and Laki suggest a common origin of the magma for this volcanic system, over at least the last eight centuries [START_REF] Bindeman | Time constraints on the origin of large volume basalts derived from O-isotope and trace element mineral zoning and U-series disequilibria in the Laki and Grímsvötn volcanic system[END_REF]. In terms of ground deformation, the present GPS array around Grímsvötn comprises a central station (GFUM) that clearly detects deformation of volcanic origin, up to fourfold greater than the expected tectonic or isostatic component, and it might be due to both shallow and deep pressure changes. Apart from this central station, due to the Vatnajökull ice cap, there is a gap of measurements between 3 km and 50 km from the caldera. The horizontal and vertical displacements that would be produced by the two sill-shaped reservoirs model are shown in Figure 10 with H d = 10 km and in supporting information Figure S3 with H d = 30 km. At a horizontal distance of 50 km from the sources, for a deep reservoir at 30 km, the horizontal displacement rate can reach 1.2 mm/year. This rate is about 25% of the horizontal GIA value we use, but it does not interfere with our results due to the small impact of GIA on the horizontal tectonic rate correction at GFUM (see supporting information Figure S3). However, such possible magmatic effects should be taken into account for the interpretation of horizontal displacements around Vatnajökull ice cap. The only signal when probing for a deep reservoir that we might expect to measure would be the exponential short-term deflation induced by the transfer of magma to the shallow reservoir following the eruption. However, Figure 10 shows that for a 10 km deep reservoir, this deflation might only be observed 5 km from the caldera center, with an amplitude of 1.1 cm about 9 months after the eruption. Due to the present poor spatial resolution, it is not possible to detect any deflation that might be produced by a deep source between 10 km and 35 km in depth. It would thus be interesting to install some more instruments on nunataks on Vatnajökull ice cap, to be able to confirm the presence of a deep reservoir and to improve our understanding of the dynamics of this volcanic system.

The model proposed here might appear oversimplified with regard to the rheology used and can be improved to account for viscoelastic behavior, at least around the deep reservoir, and for magma . Temporal evolution of modeled posteruptive vertical and radial displacement calculated at various lateral distances (r = 3.5, 5, 10, and 50 km) from the volcano axis when two sill-shaped reservoirs are considered. The deep reservoir is located at 10 km in depth with a radius of 2.2 km, and the shallow reservoir at 3 km depth with a radius of 2 km. The time constant is here set at 0.33 year, as inferred for the 2004 posteruptive displacement. The transient vertical deflation induced by the withdrawal from the deep reservoir can only be detected by a station located at 5 km from the volcano axis, where 11 mm subsidence is expected 9 months after the eruption (arrow).

compressibility. However, it has the advantage of providing a consistent explanation for the temporal evolution of the surface displacement recorded above Grímsvötn Volcano and of providing constraints on the lateral extension of the deep magma reservoir. This deep reservoir has to be of limited size (radius, <10 km), based on the vertical displacement measured (see Figure 9). Moreover, as shown in Figure 8, a too large ratio (a r ∕a s ) between the radius of the deep reservoir and the shallow reservoir would drastically reduce the transfer of magma to the surface.

Conclusions

In this study, we propose the first two-magma-chamber model that integrates magma flow and crustal deformation, and this is dedicated to interpret the temporal evolution of surface displacement in a volcanic area. The model consists of two connected magmatic reservoirs on the same vertical line, of either spherical or sill-like shape, embedded in an elastic half-space, with the magma assumed to be an incompressible fluid. This model enables the interpretation of the time series of the displacement measured by GPS at Grímsvötn, Iceland, since 1998, over three successive eruptive cycles. We derived a bottom magma production rate of between 0.01 and 0.11 km 3 /yr, which is in good agreement with previous observations, and we can conclude that the deep-seated magma reservoir is of limited size, with a lateral extension of <10 km. Further developments should include more realistic rheologies that can take into account magma compressibility and the viscoelastic behavior around the deep magma storage zone.

⎛ ⎜ ⎜ ⎝ u v w ⎞ ⎟ ⎟ ⎠ = ( 3M o 2G𝜋 ) ⎛ ⎜ ⎜ ⎜ ⎝ xH 2 R 5 yH 2 R 5 H 3 R 5 ⎞ ⎟ ⎟ ⎟ ⎠ , ( B1
)
where M o is the moment. Using this relation, the horizontal component of the surface displacement at distance r (see Figure 6) is given by

U R = 3M o 2𝜋G rH 2 R 5 . ( B2
)
To express M o as a function of the magma overpressure, and knowing that the volume of surface uplift (integral of the surface vertical displacement) for a pressurized sill is equal to the cavity volume change [START_REF] Lisowski | chap. Analytical volcano deformation source models Volcano Deformation[END_REF][START_REF] Fialko | Deformation due to pressurized horizontal circular crack in an elastic , with applications to volcano geodesy[END_REF], using equations (B1) and ( 7), we can deduce that

M o = 8 3 (1 -𝜈)a 3 ΔP. ( B3
)
It then follows that

U R = (1 -𝜈) G 4rH 2 a 3 ΔP 𝜋R 5 (B4)
and the vertical displacement )

U V = (1 -𝜈) G 4H 3 a 3 ΔP 𝜋R 5 . ( B5

Figure 2 .

 2 Figure 2. Velocities observed at the GPS stations from 2006 to 2011. (a) Horizontal component determined from the reference frame known as "ridge fixed. " (b) Vertical component in ITRF 2008. The locations of the Laki fissure (La.) and Grímsvötn volcano (G.) are indicated. For the horizontal component, a divergent direction from the fissure swarms (thin black lines) is clearly observed.

Figure 4 .

 4 Figure 4. GPS time series (three components) at GFUM station. (a, b) The horizontal components were calculated with the reference frame as "ridge fixed" and are corrected for the tectonic as explained in the text. (c) The vertical component was calculated on ITRF 2008. Vertical dotted lines indicate the two eruption dates. The reference time is the 2004 eruption, and the reference position is the last observation before this eruption. (d) The north component is shown as a function of the east component. The coordinate origin is the position of the last observation before the 2004 eruption.

  estimated from March 2006 (GFUM setup) to May 2011, beginning of the eruption. The observed vertical component at GFUM is poorly constrained due to seasonal effects and about 50 mm/yr. It may include the vertical component of the isostatic rebound estimated by Geirsson et al. [2010] at 22 mm/yr.

Figure 5 .

 5 Figure5. Comparison of the ground deformation observed at Grímsvötn during the last three eruptive sequences and corrected from the estimated tectonic component. (a) Directions of posteruptive displacement, up to 6 months after the eruption, and coeruptive displacement, recorded at the closest GPS station from the volcano, GRIM (replaced by continuous GPS station GFUM in 2004). The Grímsvötn ice surface map is modified from[START_REF] Gudmundsson | Eruptions in Grimsvötn, Vatnajökull, Iceland[END_REF], and the locations of the last three vents are shown (the 2004 and 2011 vents are the same). We observe the same direction for the posteruptive and coeruptive displacement, even if the eruptive vents were not located exactly in the same place. (b) Radial displacement measured at the GRIM station (after the 1998 eruption) and GFUM station (after the 2004 and 2011 eruptions). The fits for equation (20) are shown as dotted lines.

  Figure10. Temporal evolution of modeled posteruptive vertical and radial displacement calculated at various lateral distances (r = 3.5, 5, 10, and 50 km) from the volcano axis when two sill-shaped reservoirs are considered. The deep reservoir is located at 10 km in depth with a radius of 2.2 km, and the shallow reservoir at 3 km depth with a radius of 2 km. The time constant is here set at 0.33 year, as inferred for the 2004 posteruptive displacement. The transient vertical deflation induced by the withdrawal from the deep reservoir can only be detected by a station located at 5 km from the volcano axis, where 11 mm subsidence is expected 9 months after the eruption (arrow).

Table 1 .

 1 Results for the Modeled Tectonic Trend and the Residual Volcanic Trend at the GFUM GPS Station, Estimated From the Post-2004 Eruption

	GFUM Station	North Component East Component	Amplitude
	GPS data	-38.7 mm/yr	29.3 mm/yr	48.5 mm/yr
	Tectonic component estimated by model	-2.7 mm/yr	7.5 mm/yr	8.0 mm/yr
	Residual trend	-36.0 mm/yr	21.8 mm/yr	42.1 mm/yr

Research: Solid Earth

  Illustration of the model geometry with two magmatic reservoirs. The first shallow sill-like reservoir at depth H s and radius a s is fed by magma that flows through a cylindrical conduit of radius a c and length H c . The second reservoir is deeper, of radius a d , and located at depth H d . Both of these magma chambers are located under the volcano at the respective distance R s and R d from a GPS station, named S. The radial distance at the surface between the volcano axis and the GPS station (S) is denoted as r.

	10.1002/2013JB010569

Figure 6.

Research: Solid Earth

	10.1002/2013JB010569

Figure 7. Evolution of overpressures ΔP s (t) and ΔP d (t)

Journal of Geophysical Research: Solid Earth

  UV ∞ , ΔV, and Q ∞ . The effects of the error on the shallow reservoir depth H s are less straightforward to calculate, and the induced errors on ΔV and Q ∞ can be quite large and will depend on H s . A change from 2.5 km to 3.5 km in H s will increase UV ∞ by up to 35%, due to the deepening of the ground deformation source. A shallow reservoir deeper than 4 km is not possible at Grímsvötn, due to the low ratio of the vertical/ horizontal UV ∞ / UR∞ displacement trends, which is lower than 1.2. It is also of note that the Solutions of the model producing the observed rate of horizontal displacement UR∞ at GFUM for different depths H d (in km) and reservoir radii ratio (a d ∕a s ) (see Figure

	10.1002/2013JB010569

Figure 9.

©2014. American Geophysical Union. All Rights Reserved.

Acknowledgments

We gratefully acknowledge Halldór Geirsson, Sigrún Hreinsdóttir, Richard A. Bennett, Thóra Árnadóttir, Benedikt G. Ófeigsson, and Thorsteinn Jónsson, who supervised and led the installation and operation of the high-rate GPS stations used in this study. This work was supported by the IPEV, Jules Verne grant, and Université de Savoie. Also, we thank Freysteinn Sigmundsson for fruitful discussions, the Associate Editor, and two anonymous referees for their constructive comments and suggestions.

Appendix A: Resolution of the Differential Equation in 𝚫P s and 𝚫P d

It follows from equations ( 4), (6), and (8) that

In the same way, using equations ( 4), (5), and (7), we get

We thus obtain

Combining equations ( A3) and (A4), we obtain

By integrating equation (A5), we get

where ΔP d (0) and ΔP s (0) are the initial overpressure in the deep and shallow reservoirs, respectively.

Appendix B: Calculation of Horizontal Displacement at the Surface Induced by a Sill-Like Magma Reservoir

The three components of the displacement vector induced at the surface at point M(x, y, 0) by a sill-like magma body embedded in an elastic half-space is given by [START_REF] Lisowski | chap. Analytical volcano deformation source models Volcano Deformation[END_REF]: